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SUMMARY 


Analytic models are developed for computing the periodic sound pres- 
sures of subsonic fans and compressors in an infinite, hardwall annular 
duct with uniform flow. The basic sound-generating mechanism is the scat- 
tering into sound waves of velocity disturbances appearing to the rotor or 
stator blades as a series of harmonic gusts. The models include the sources 
of velocity disturbances arising from the component interactions studied 
by Kemp and Sears, and non-component-induced inlet distortions at a rotor. 
The primary result of the models is the computation of the periodic sound 
pressure mode amplitudes which can be used as inputs to a duct acoustic 
program for computing the propagation through a finite, multisectioned 
duct and the radiation to the far field outside the duct. Computer sub- 
programs for the mode amplitude calculations are presented. 
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1.0 INTRODUCTION 
1.1 Background 

Gutin (ref. 1) in 1936 published the first analytic model for the 
tone noise generated by rotating machinery. His idea was to replace a 
rotating propeller by a disc of applied acoustic dipoles, the strengths 
of which were determined by the steady air forces acting on the propeller 
blades and the frequencies of which were determined by the propeller shaft 
frequency and the number of blades. Not until 1963, when Van de Vooren 
and Zandbergen (ref. 2) published their work, did the problem receive an 
essentially different formulation: The rotating and translating blades 

were replaced by dipoles and monopoles , the strengths of which were to 
be determined by solving a linear boundary value problem. They did not, 
however, solve this problem, but returned, instead, to the momentum con- 
siderations of Gutin for the strengths of these moving sources. 

The advent of the turbofan engine gave a new impetus to the construc- 
tion of analytic models for the generation of tone noise by rotating ma- 
chinery. In 1961 Slutsky, et al. (ref. 3) developed analytic models for 
the generation of tone noise by a turbofan stage assuming for the sources 
of sound the steady, rotating vorticity of the rotor blades and the 
unsteady vorticity of the rotor and stator blades induced by the two- 
dimensional cascade aerodynamic interactions studied by Kemp and Sears 
(refs. 4 and 5). These noise models were for a ducted fan stage and 
included, in principle, the transmission and reflection at the termina- 
tion of the duct. Tyler and Sofrin (ref. 6) reported upon similar, 
althrough less detailed, work, e.g., no explicit source modeling, in 
1962, and a year later Hetherington (ref. 7) reported upon a source 
modeling scheme essentially the same as that used by Slutsky (whose 
work did not receive general notice until much later) . Griffith (1964, 
ref. 8) added randomness to the phasing considerations in the blade row 
interaction tone noise generation processes and Hulse, et al. (ref. 9) 
added similar, although less detailed, considerations to the Gutin-type 
of free-rotor tone noise mechanism. At this time Lighthill's (ref. 10) 
theory of aerodynamic noise, and Curie's (ref. 11) extended version, 
began to be applied to the problem of the analytical modeling of the 



tone noise generation by fans and compressors, particularly by Lowson 
(refs. 12 and 13. Ffovcs Williams and Hawkins (ref. lb) advanced an 
argument based on the Lighthill theory for a tone noise generating mech- 
anism for a fan or compressor rotor that was neither a surface dipole or 
monopole, but a fluid volume quadrupole . They did not, however, develop 
the idea into a working analytic model. Wright (ref. 15) and Horse and 
Ingard (ref. 16) developed working models for tone noise generation by 
free rotors in steady nonuniform flow using time- dependent Gutin-type 
dipoles, and Lowson (ref. 17) added to these models the element of non- 
steadiness in the nonuniform flow at the rotor (see Griffiths and Hulse) . 
Pfenninger (ref. 18) gave physical insight to the problem of tone noise 
generation by turbofans in non-steady, non-uniform inflow conditions. 

By employing an interesting combination of Griffith's and Pfenninger 's 
ideas, Hanson (ref. 19) has recently advanced a model (without explicit 
source modeling) for computing the complete spectrum of sound radiated 
by a free rotor with turbulent inflow, and a similar blade row inter- 
action model (refs. 19 and 20, with explicit source modeling). Detailed 
source modeling for the tone noise generation by a ducted rotor in uni- 
form flow is being investigated by Lordi (ref. 21) while Drischler 
(ref. 22) has computed pressures in the duct of a ducted propeller using 
a Gutin-type dipole model for the blades. Schaut (ref. 23), in using the 
Morse and Ingard radiation model for a free rotor, modeled the blade 
pressures after the design profiles. 

1.2 Present Work 

It was the purpose of the present work to develop an analytical 
procedure and a set of computer subprograms for computing the sound pres- 
sure at the harmonics of the blade passing frequency which is radiated 
from an axial flow compressor or fan stage in an infinite, hardwalled 
annular duct. The problem is presented in the form of a boundary value 
problem to which approximate solutions are adapted from the published 
literature on the aerodynamics of thin airfoils. The radiated pressure 
in the duct is most conveniently expressed in terms of an eigenfunction 



expansion, with the coefficients of the expansion determined by the 
pressure on the fan rotor and stator blades. By the familiar process 
of linear analysis, these coefficients then can be used as inputs to 
problems involving a large variety of ducts. Zorumski (ref. 24) has 
recently proposed a very extensive and systematic scheme based on mode 
matching techniques, which utilizes these infinite-duct coefficients 
as the elements of an input, or known, vector in an inhomogeneous matrix 
equation wherein the unknown vector is the set of coefficients for the 
pressure expansion in the eigenfunctions appropriate to the section of 
duct contiguous with the inlet or exit of a finite, multisectioned duct, 
which, when known, allows the computation of the far-field pressure out- 
side the duct (see, e.g., Slutsky [ref. 3], Lansing [ref. 25], and Clark 
[ref. 26]). The inverse of the matrix multiplying the unknown vector 
represents the multiple transformations which the pressure expansion 
coefficients undergo between the duct section containing the compressor 
or fan, and the section contiguous with either the inlet or the exit. 

It was, in part, to supply computer subprograms for computing input 
vectors to this matrix equation that the present work was performed. 
Examples of this matrix equation in the present context can be found 
in Zorumski' s report and Lansing and Zorumski (ref. 27). 

The models considered include the blade row interactions of Kemp and 
Sears and the rotor alone in steady and quasi-steady nonuniform inflow. 
The computer subprograms have been coded in standard FORTRAN for a CDC 
6600 scientific computer. They have been designed to compute the non- 
dimensionalized, infinite duct coefficients, or mode amplitudes, for 
the eigenfunction expansion of the radiated pressure for each of the 
models. Since these are subprograms and their outputs are intermediate 
results in the computation of the radiated pressure, a user must supply 
a main program for computing pressures from the expansion coefficients. 



2.0 ACOUSTIC THEORY 


This analysis is based upon the linearized theory of compressible 
fluids, thus restricting the applicability of the results to lightly 
loaded blades and to subsonic subcritical relative flow velocities. 

Some of the tone noise characteristics of modern transonic fans, e.g., 
the "buzz saw" phenomenon, are, therefore, excluded from consideration. 

It is assumed that under these conditions the most efficient discrete 
frequency sound-generating mechanism is the scattering into outgoing 
pressure waves of flow perturbations which appear to the turbomachinery 
blades as unsteady velocity disturbances. Since a perturbation-free 
flow is difficult to achieve under the best of conditions, this mechanism, 
it seems, will persist into the flow regimes in which, presently, dynamic 
calculations cannot be made. Under linear conditions, the strengths of 
the outgoing pressure waves are proportional to the strengths of the 
velocity disturbances at the fan rotor and/or stator faces. The analysis 
consists of determining the transfer function for general velocity dis- 
turbances and then identifying and describing the velocity disturbances 
which occur in the operation of ducted turbomachines. Although little 
work has been done in the determination of the transfer function, even 
less has been done in analyzing and describing these velocity disturbances. 

To obtain the transfer function requires doing two things: solving 

the convective wave equation with a dipole singularity as an inhomo- 
geneous term, and then solving for the induced dipole distribution that 
is required to cancel out the incident velocity perturbations that are 
normal to the rotor and/or stator blades, which are taken to be thin 
airfoils. These are discussed in this section, leaving to the next 
section the details of the tone noise models that result when the 
incident velocities are specified. 


2.1 The Acoustic Propagator 

The geometry of the problem is depicted in figure 1. The mathemat- 
ical conventions and the nondimensionalized form of the linearized fluid 



Annulus formed by 



cylinders. 





CTION A-A 




equations are presented in appendix A. The governing equation for the 
linear pressure variations in a fluid with uniform mean flow velocity, 

M, parallel to the positive z-axis of the coordinate system, that results 
from the presence of a dipole of strength -F(t), orientation 


~ | o,e^,ez | and position r s = | P s* ^ 3 , z sj 
-□ c P(r,t) = F(t )e4 6 (r - r ) 


is 


( 2 . 1 . 1 ) 


where □ is the convective D* Alembert i an, 
c 


^ « £ ) 2 - 2 
~\at 3z ] p 3p 


( 2 . 1 . 2 ) 


3p 


_1 __3 

p2 


3 3z 


and V is the gradient with respect to the s-subscripted coordinates, 
s 


If a solution to the simpler equation 




□ T(r,r ,t-t ) = 5 (r-r )6(t-t ) 


(2.1.3) 


is found, then the solution to equation (2.1.1) is 

CO 

p(?,t)= / F(t o )e4 o r(?,? o , t-t o )A = + dt Q 
-°° /os 


(2.1.4) 


In this report T will be referred to as the acoustic "propagator," 
a descriptive word often used in similar contexts in elementary particle 
physics, e.g., Feynman (ref. 28). It satisfies causality but, in this 



case, not reciprocity, so that it is not a Green's function in the com- 
monly accepted sense, e.g., Morse and Feshbach, chapter 7 (ref. 29). 

For equation (2.1.4) to hold, T must satisfy the same boundary conditions 
as does p. These conditions are that the normal gradient of p at the 
inner and outer duct walls should vanish. 


3P 

8p 


o at p=l, p=r). 


(2.1.5) 


and that only outgoing waves should exist. The problem is simplified 
by defining the temporal Fourier transform of T by 


A 

r 


oo 

(r,r o ,w) = / r (r,r o ,t)e dr 

— oo 


( 2 . 1 . 6 ) 


and the "generalized Prandtl-Glauert" transformation by 


A _ ^ 
r (r,r o ,w) 


.. Mo)(z-z 0 ) 

„ " 1-M2 
e 

V l-M 2 



(2.1.7) 


where 



( 2 . 1 . 8 ) 


and 


r* = |p,<),z*| 


with 



(2.1.9) 


7 



Then the equation y must satisfy is 


T * 2 a 
V y + u y 




( 2 . 1 . 10 ) 


with the boundary conditions that 


/\ 

ll 

3p 


o at p-i, p»n, 


( 2 . 1 . 11 ) 


and for | z-z | a ** constant, 

s 




( 2 . 1 . 12 ) 


This last condition is the statement that if the fluid is considered 
to have a small absorptivity, equivalent to letting o>* have a small 
imaginary part, then the pressure decays exponentially away along the 
axis. However, when the absorptivity is neglected, Im(io*) * 0, then at 
least a part of the pressure does not decay for | z-z I-* 00 (neglecting, 
of course, all other absorption mechanisms) . It is convenient to solve 

a 

for y as the sum of a particular integral to the inhomogeneous equation 
neglecting the presence of the duct walls, with the general solution to 
the homogeneous equation, and then adjust the undetermined parameters in 
the general solution to satisfy the hardwall duct boundary conditions. 

This is done in appendix B, from which equation (B42) is found: 


QO 00 

y - (uo)fl (y ,P_)d (z* - z*,oj*) (2.1.13) 

2**-' m mn b'Vo' mn o’ 

m=-°° n= 0 


8 



with the definitions as given in appendix B. This result could be re- 
ferred to appropriately as the Prandtl-Glauert space Green's function. 
The acoustic propagator is, then, from equation (2.1.7): 


I**—* £■— m mn m ran o 

m » n “o 


* ®( z " z o»“) 
mn 


(2.1.14) 


with 


..-Min / x 
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J mn - 


/l-M 2 


j / k v 
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mn o 


.Mm 


_ ± ^ (z _ ) . 


$ 
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= e 


1-M 


1-M 

2ig 


z-z 


mn 


(2.1.15) 


where 


0 


mn 


-VI 


2 2 
u) -(1-M )p 


2 

ran 


(2.1.16) 


and 8 has the same properties as a function of u that s has a func- 
mn mn 

tion of u* (see appendix B) , but with different branch points. 


An alternative way of expressing D is 

mn 


D 


mn 


D + 
mn 


+ D_ 


mn 


(2.1.17) 
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with 


( 2 . 1 . 18 ) 


«+ i^mn (z-z ) 

D~ = e(±(z-z )) — 

mn o " 2i8 

mn 


e(X) = 


(l,X>o 

[o,X<o 


and 


-a)M±f3 

+ mn 


mn 1-M 2 


This makes explicit the distinction between downstream (+) and 
(-) propagation. 


Finally, the time-dependent propagator is given by 

T (r ,r 0 ,t) = P)^ m (P^P °> D mn (z-z 0 ,x) 

On «— * m mn m mn mn 


m=-°° 


n = -° 


with 


D 

mn 


1 

2tt 


a> 

r d 

J mn 


(z-z 0 ,a))e'^ i)T dci) 


= D+ + 
mn 


D- 

mn 


and 


D± 


mn 


z-z 0 

e(x± 

1±M 




(z-z 0 -Mx) 


2 


) 


( 2 . 1 . 19 ) 

( 2 . 1 . 20 ) 

upstream 

( 2 . 1 . 21 ) 

( 2 . 1 . 22 ) 

( 2 . 1 . 23 ) 
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This result follows from the discussion in Morse and Feshbach, chap- 
ter 2, (ref. 29), on the Klein-Gordon equation and the elastically braced 
string. This reference also is useful for understanding how the chan- 
neling by the duct walls results in the fluid being a dispersive medium 
for acoustic propagation. Otherwise, this result for D can be found 
by converting the Fourier inversion to a Laplace inversion for causal 
functions, and looking up the result in a table of Laplace transforms 
(e.g., equation [29.3.92], p. 1027 of Abramowitz and Stegun, ref. 30). 


2.2 Radiation from Dipoles in Ducts 

■V 

The solution to equation (2.1.1) for a stationary dipole (r 
independent of time) is, from equations (2.1.4) and (2.1.14), 


1 00 A 

p(r , t) f p(r,u>)e awt du) 

2ir J 


where 


-> -± 


p(r ,oi) = F (m) e* V T (r , r ,u>) 



£ [ e(z-z s )A ^ n (w)e iK mn Z 
n=o 


+ e(z -z)A-(a))ei K ^n z 1 d (y p ) 
s m n m mn 


( 2 . 2 . 1 ) 


with 


A± 

mn 


(u) 


-F(u») 


g^e +K± e 
H s ji an ; 


23 


mn 


\ n . N -x(m d> + K± z ) 

(ft (y p )e Y s mn s 
1 m mn s 


( 2 . 2 . 2 ) 
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These are the downstream or upstream mode amplitudes for a station- 
ary dipole source. Since F(t) is a real function, it follows from the 

definition of 3 (see the discussion of s in appendix B) that A - (u) 
mn mn run 

has the symmetry property 


A -mn (_a>) = (A Sn (a)))+ * 

-mn mn 

allowing negative frequency mode amplitudes to be known in terms of 
positive frequency mode amplitudes. This symmetry guarantees that 
p(r,t) is real, and 


00 p 

p (r,t) =2^-/2 REAL [p(?,m)i iait dm 


(2.2.4) 


When the same dipole is rotating at the constant rotational fre- 
quency £2 while maintaining its orientation with respect to the local 

cylindrical coordinate triad, its position vector, r , becomes time- 

s 

dependent, 


r 

s 



with 


(t ) = fit+(j)-27r ,£ = —<*>,•• «o, •** , + 00 (2.2.5) 

where i l is the revolution counter and <j> is the angle coordinate of the 
dipole for £ = 0 and t = 0. Then, 


oo 

F(t) — ► F ( <f> ( t) ,t) = F(t)^p6(<J>-ftt + $+2tt£) 

£ = — 00 


( 2 . 2 . 6 ) 



and Che pressure field of this rotating dipole is given by 


P(r ,t) 


® 2 * 


j J F (W' t o) 


e*V T(r,r ,t-t ) 
o o o 
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dd> d t 
p o o 
s 
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= 4— ^ y 61 (y p) & ( u P) e im< ^ 

2tt / — > / ^ m mn m v mn o J 
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03 


00 2ir 





O 


(♦o^V* + 2ir£)F(to) 


im 

P c 


D (z-z t) + 
mn & 


9D (z-z ,t) 
mn s 

9z 


imd> 

o 

e 

z 

— 


d<t> dt 
o o 


(2.2.7) 


The evaluation of this expression is performed in appendix C. The 
result is, for the pressure spectral density downstream and upstream, 
(±) , respectively, 


~±-y 

p(r,w) 


2ir 


X—i mn m mn 


i(md» + K± z) 
mn 


m=-°= n=o 


( 2 . 2 . 8 ) 
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with 


A± (u) = -F(cj-mQ) 
mn 


m 

p 6 <J) +K± e 
s mn z 

2B 

mn 


a/ \ i(m<j) + K± z ) 

(R (p p )e y mn s'' 
m mn s 


(2.2.9) 


The same symmetry property prevails, as it must, but there are two 
differences from the mode amplitudes for a stationary dipole: <£ is the 

angle coordinate for the dipole location in the reference frame rotating 
with the dipole, and the dipole spectral density is shifted differently 
for each spinning mode (m is the spinning mode index). 

The pressure field of a stationary surface distribution of dipoles is 
the sum of the individual dipoles associated with each infinitesimal 
area of the surface. Letting 


F(«o)-F(p st 5(+ s ,* s ),w) 


( 2 . 2 . 10 ) 


be the dipole surface density spectral density, where it is assumed that 
the surface extends from the inner to the outer duct wall and for each 
radius the lateral extent of the surface projection onto a duct cross 
section is sufficiently small so that a straight line can replace the 
arc (see fig. 2). The surface can then be taken to be made up of dif- 
ferent straight line segments having projections on <p and z, but not 
on p which are strung together from inner to outer duct radius, i.e., 
a fan blade whose surface is, at each span segment, approximated by its 
chord line. If e = jo, e , is the unit vector normal to the surface 

for a given radius, then the local rectangular coordinates are £ , f , 
where (see fig. 2) 


£ = -p <P' e 
s s z 


+ z 


e 

s <p 


C 


p <i 1 e + z 1 e 
S T S (j) 8 z 


( 2 . 2 . 11 ) 
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FIGURE 2. - DIPOLE SURFACE GEOMETRY 



and 


with 


z' = ?e + |e, 

s z <J> 

p d>' = Ce,-Ce 

H s Y s <p z 


*s = 4, 


M.C. 


+ 


z 
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z 


M.C. 


+ z' 

s 


( 2 . 2 . 12 ) 


(2.2.13) 


where C * Z M C are t * ie m ^dchord coordinates at the radial station, 

p g . The phase factor in the expression for A“ n in equation (2.2.2) 
becomes 


-i (m<)> + K± z 

s ran 


) -i (m4> +K± z M _) 

51 ' M.C, mn M.C. 

= e 


-i ( K± e e )£ 
mn 4> P„ z 7 ^ 
e *» 


2.2.14) 


when evaluated at f =0, i.e., on the chord line. Hence, the mode ampli- 
tude to be used in the expansion for the pressure radiated by this dis- 
tribution of dipoles is 


Ajt (<»>) = 

mn 


— i (rntf).. + K± z 

e ' M.C . mn 


28. 




1 

( r e m 

J ♦T 


+ e K+ 
z mn 


<R (p P ) 

m ran s 


* j F(p s ,C,u)e" i ( K &nVPs e z ] 5 d£dp g (2 215) 

-c/2 

w 7 ith c the lateral extent of the surface at the radius, p g , i.e., the 
chord of the equivalent blade segment . When there are N such surfaces 
equispaced about the annulus with the same z r but with 

iM ■ L* • 


*M.C.-*J 



(2.2.16) 
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and 


F 




then 


A± (w) 
mn 


N-l 


*Z) A Sn3 ( “ > 


-iK± z 
mn 


28 


mn 


-/ 
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♦p. 


e K± 
z mn j 


iR 
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^ v mn P s^ 


(2.2.17) 




e’ 


d5,dp s 


with C ’ » S/c/2 and 


+ 

1C «* 

mn 



(2.2.18) 


Similarly, for a rotating set of surfaces 
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mn 


e ^mnSl.C. 
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1 N-l . 2-rrj 

* f f { £ e‘ N F, (p s , 5 ’,»-4 
-1 j=o 

*e 1<mn5 dC.'dps (2.2.iy) 

2.3 The Dipole Surface-Density Function 

The dipole surf ace -density functions are the pressure difference 
functions of thin airfoil theory and must be determined from the condition 
that the normal velocities at the airfoil surfaces vanish. When the blades 
comprising a fan or compressor blade row are unstaggered, flat plates, a 
precise formulation in terms of a set of integral equations can be made 
(see appendix D) . A similar formulation for a blade row which performs 
a net turning of the mean flow is as yet unavailable (see Atassi and 
Goldstein, ref. 31, for a linear formulation for the unsteady lift with 
finite camber and angle of attack of a single two-dimensional airfoil) . 

It was not, however, the intention of this work to develop this formula- 
tion and seek solutions within it. The intention was to recognize the 
basic requirements for such a formulation and seek means of approximately 
satisfying these requirements. The basic requirement is that the incident 
velocity perturbations normal to the blade surfaces should be cancelled 
by the induced velocities of the ducted blade row. The approximations 
to the satisfaction of this requirement are determined by the available 
aerodynamic calculations and the form of the result of these calculations. 
It is desired that the result be in the form of known mathematical func- 
tions or numerical routines the evaluation of which is not much more 
cumbersome than the evaluation of known functions. A further require- 
ment is that the approximations should be consistently employed. The 
approximation used in this work is that the incident normal velocities 
are cancelled by the velocities induced in an incompressible fluid 
neglecting three-dimensional effects, the presence of the other blades 
in the row, and the duct walls (see appendix D) . In particular, the 
calculations of Sears (ref. 32), Kemp (ref. 33), Horlock (ref. 34), 

Naumann and Yeh (ref. 35), and Filotas (ref. 36) are used. The mean 
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flow conditions are assumed to be those through the staggered row (see 
fig. 3). The general form for the dipole surface-density function is, 
then, 


(P„»£',w) = 


dCr 


do 


u ^ P s. Z s,^ - H (p c 


s. 


(2.3.1) 


mean nondimens ionalized velocity through the row, 

spectral density of the normal velocity perturbation 
th 

incident on the j blade 

nondimens ional dipole surface-density response function 
for convected, harmonic gusts 

Of the models to be considered, only the potential field blade row 

interaction model does not make use of the assumption that the incident 

velocities can be taken to be "frozen-convected gusts." Under this 

assumption, u. (p , z , u) is taken to be independent of £' , i.e., its 
J s s 

variation over the chord of the airfoil is negligibly small and z - z , 

S M • u • 

or some other appropriate value such as the axial coordinate of the 
point on the airfoil one-quarter chord from the leading edge (see Sears, 
ref. 32). The item dC L/ da i s the slope of the steady lift versus angle 
of attack curve, which for a flat, two-dimensional airfoil is 2tt . The 
use of dC L/ da instead of 2-rr arose from the quasi-steady lift consider- 
ations: if low-frequency sound pressure is proportional to the section 

unsteady lift, and 


where: 

= 

u . = 

J 



r r 


1 

IT 



-1 


C’.u) 


d£' 


1, 


UKO 
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then 


<$L 


1 2 dC L 

2 cm m ST { » 


with 



where : 
6L 
6a 


section lift variation 
angle of attack variation 



3.0 MODELS FOR TONE NOISE GENERATION 


The four models discussed in this section are referred to as the 
viscous wake interaction model, the potential field interaction model, 
the rotor in steady distortion model, and the rotor in nonsteady dis- 
tortion model (or the "single eddy" model) . A model consists of de- 
fining the u in equation (2.3.1) along with appropriate H. Lifting- 

J -4- 

line models result when the < in the chordwise integration is taken 
to be zero. The potential field interaction model is a lifting-line 
model, while the other models include the lifting-line assumption as 
an option. Assuming a lifting line is equivalent to assuming the blade 
is a chordwise compact acoustic source. It is felt, in turn, that a 

first order model for a chordwise noncompact airfoil is achieved by net 

+• 

setting equal to zero. Higher order estimates will include com- 

mn 

pressibility in the dipole surface-density function, which means 
including the effects of the other blades in the row and the duct walls 
as well as the propagation of pressure waves from one part of the air- 
foil to another. 

The specific results of each model will be an analytical procedure 
for computing the mode amplitudes of the pressure radiated by a single 
component of a fan stage, i.e., equation (2.2.17) or equation (2.2.19), 
with a specification of the chordwise integral 


c_ 

2 





F. <P S ,€'.«> 


-IK" S' 

, 11111 dg' 


(3.0.1) 


or, in the lifting-line models. 


N-l . 2-tt j 

-lm ■ v a 

N T . , 

e Lj (P a ,w) 


E 


1=0 


(3.0.2) 



where 


Lj (P s ,m ) 




(P s .5' ,«*>) 


&v 


(3.0.3) 


3.1 The Viscous Wake Interaction Model 

In an ideal axial flow fan or compressor stage, the inlet flow and 
exit flow are axial. To achieve this requires at least two components — 
a rotor and a stator. Viscous wake interactions occur when a rotor cuts 
through the viscous wakes from the blades of an upstream stator and when 
the viscous wakes from an upstream rotor sweep by the blades of a stator. 
These possibilities are given a two-dimensional cascade representation 
in figures 4 and 5. Consider first the inlet stator-rotor configuration 
(fig. 4). Since the wake defect will be known in the primed coordinate 
system attached to the inlet stator blades, it is necessary to determine 
the transformation from this coordinate system to the system attached to 
the rotor blades (the "bar-primed" coordinates of fig. 4). This trans- 
formation is (subscript 1 is for upstream row and subscript 2 is for 
downstream row) 

Z' = V COS 0 + Y* SIN 6 + d COS ip + UtSIN 

(3.1.1) 

Y* = Z* SIN B + 7' COS B - d SIN ip + UtCOS ifi 

t It t It 

Then the coordinates for the £ wake in the j rotor blade 
coordinate system are 

Z' = Z' . COS 0 + Y' . SIN B + d COS \p + UtSIN ip + ( j b^-iUm) SIN ip 
& 3 J 3 

Y' . =-Z' . SIN 0 + 7' . COS 0 - d SIN <p + UtCOS + ( jb - Zb.) COS <p 
J J 2 1 

(3.1.2) 

where all the symbols correspond to those used in figure 4 and • 
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FIGURE 5. - CASCADE REPRESENTATION OF ROTOR OUTLET STATOR COMBINATION 
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u = 

rotor blade nondimensionalized speed equal to M^p 

M = 

T 

rotor blade tip nondimensionalized speed 
g = relative stagger angle, y + 

Y = rotor blade row stagger angle, 
ij, *= exit angle of the stator row. 

If is the number of blades in the upstream row, the stator row, and 
N 2 is the number of blades in the downstream row, the rotor row, then 


b^ = 2-irp/N^ 


and 


b 2 = 2trp/N 2 


(3.1.3) 


. th 


Representing the wake defect velocity behind the Z 
in the j 1 "* 1 rotor blade coordinates by W and assuming it can be 
developed in a double spatial Fourier integral, then 


stator blade 



CO 

AA 


dA dK=. 

2tt 2 it W(K,X)e 


i (KY' + AZ’ ) 






(3.1.4) 


and the velocity due to all the stator blades is 

OO 

w - w„ . 

3 L-j £j 

£=-00 


(3.1.5) 


00 


dK = ™ i(KY' + AZ’..) 

dx dK W(K>X) yy 13 £j 

£=-co 


2tt 2ir 
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(3.1.5) cont. 


so 
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iK [-ZV SINg + Y'^COSg - dSINtp + UtCOSt//| 




iA f Z * .COSS + Y' . SINg + dCOS^ + UtSIN^I 

*e L J J J 


i KCOSip + A Silty 3b 


j| 31 
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J£=-oo 
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-Ub, KCOSij/ + ASItty 


From the Poisson sum rule 

-ifi.b 1 J^KCOSi^ + XSINxpJ 




il=-oo 


CO 

= b^COS^ X] 6 [ K "(b^COS^~ ATAN ^ 

n=— oo 


(3.1.6) 


.th 


that the velocity at the surface of the j rotor blade is 

dTAlty 


W. ( Y* . = 0)= - 
3 3 b 


COSi|/ 5 ~J 

n = n 


b 2 

2rinj — -2rin , 
b, b 

e i e 1 


n=-co 


(3.1.7) 



(3.1.7) cont. 


o-.,- n (z'.SINg - UtCOSijj ) 
*e b-jCOSiJ; 3 
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< X ) 



2m 

bjCOSiJ) 


ATAN^ , A| 


i.\ 
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d + Z ' . COSy 
3 

COSiJ; 


If the stator wakes do not decay too rapidly along Z' , then 
W(K,A) will be concentrated in a narrow bandwidth about A = 0. Then, 
if ij> is not close to tt / 2 the tan \Ji term in the K-dependence of W can be 
dropped, giving for the A integral. 


- ATAN^e 1 


d + Z' . COSy 

A J 

COS^ 


~ f dA TT / 2irn A i*Z' 

~J 2i w fccos? ’V 

-oo ' L 


(3.1.8) 


where Z 1 is the same for all blades and is the distance from the stator 
midchord to the axial plane cutting through the rotor blade quarter 
chord (see the discussion after equation ^2.3.1]) as measured along the 
mean streamlines from the stator row (see fig. 4). 



Then, when and b 2 from equation (3.1.3) are used. 
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2tt pCOS^ W \ PC 0S ^’ 
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(3.1.9) 


-inN 

*e 
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1 pCOSij; 
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with 


UCOSg 

SIN^ 


(3.1.10) 


This velocity perturbation on the 


. th 
3 


rotor blade is thus a sum of 


f rozen-convected harmonic gusts, with u. given by 


00 



— oo 


(3.1.11) 


oo 



n=-oo 


N 

_i SINS 

p COSt|; W 
s 



e-inNi 


dTAN^ 


6 (w + nN^) 



This fora for justifies the use of the pressure difference 
function induced by f rozen-convected harmonic gusts given in 
reference 35 (see equation [26 lof this reference) for the nondimen- 
sional dipole surface-density response function (with £'= - cos 9) 


H (p ,S',a))= S(v) COT -| - COT$[aCOT f e 1VCOS0 
s 2 L 2 


(3.1.12) 


{■ 


ivcosO 


+f ^ F(v) COT j + 2 SING e 


v H lk ' lj k (v) SINk6 }] 


k=l 


with the reduced frequency defined by 


0)^2 

V “ V 2 


(3.1.13) 


and S (v) is the Sears function 


S(v) = - 


■p <-H o (2) (v) + iH x (2) (v)) 


(3.1.14) 


and F (v) is defined by 


r 

F(v) = T(v) [j*<v) - J 


- J(v) + 


J L (v) 


(3.1.15) 


with 


J(v) = J q (v) + iJ^(v) 


(3.1.16) 
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and 


(3.1.17) 


T(v) 


(v) + i H<2 )(v) 
-H^( v )+ i (v) 


( 2 ) ( 2 ) 

and Hq , are the Hankel functions of the second kind. The 

chordwise integration gives the acoustic response function of the air- 
foil span section 
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(3.1.18) 


The section lift response is from (equation [27] of ref. 35): 



(P s ,5'»“) 
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jaF a (v) + f F f (v)J 


(3.1.19) 


30 



with 


F (v) = J(v) 
a 


(3.1.20) 


ana 


J, (v) 

F f (v) = F(v) + 4 — - — 


(3.1.21) 


In these formulas, a is the small angle of attack (in radians) and 
f is the ratio of maximum camber of a slightly cambered airfoil over the 
half-chord length (see fig. 6). 

Combining these results with equation (2.3.1), the chordwise inte- 
gration in equation (2.2.19) becomes (with - M^) 
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(3.1.22) 


(See appendix E for a discussion of the d-dependent phase) . 




The sum on j can now be done: 
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(3.1.23) 


and gives new form for the right-hand side of equation (3.1.22), where 
n + ( to avoid confusion with the mode amplitude n index. 
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2M 

where V is approximated by M 2 , and g by the sum, Y + (j (see fig. 3). 



There remains to be defined the wake defect in terms of the primed 
coordinates of a typical stator blade, W(Y', Z') , in order to get 
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e- iKY ’dY' 


(3.1.26) 


The choice is made to use the empirical formulas of Silverstein, 
et al - (ref . 37) : 


W(Y',Z') = W q (Z') COS 2 



(3.1.27) 
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drag coefficient of the stator blade at P s> 


exit velocity from the stator row at p 


(3.1.28) 
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c = stator blade chord at P . 
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These formulas give for 
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(3.1.32) 


from K = 1 N,/p cos if> 
l s 


This completes the specification of the mode amplitude calculation 
for the viscous wake interaction model with a rotor downstream of an 
inlet stator. The final formula is 
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(3.1.35) 


(3.1.36) 
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The computer subprogram AAAAA (see section 3.1.1, volume II) 
computes the factor multiplying the frequency spike, 

6 ( to - c Ngl^,) , for w in the propagating region. 


With subscript 1 referring to the upstream rotor and subscript 2 
to the downstream stator , the same development for the viscous wake 
interaction model with a rotor and an outlet stator can be made, and 
gives for the mode amplitude 
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A^is given by equation (3.1.35), by equations (3.1.36) and (3.1.37), 
with >!->a , and G is given by equation (3.1.38) with 


V £ - “ V o = " 0N 1 




M. 


2M 


(3.1. 42) 


The computer subprogram AAAAA computes the factor multiplying the 
frequency spike for w in the propagating region. 


3.2 The Potential Field Interaction Model 

In the close vicinity of a blade row, the potential flow field is 
not uniform. This is due to the noncontinuous distribution of the aero- 
dynamic forces that are exerted by the blades onto the medium. If an 
airfoil moves through this nonuniform flow field, it will experience 
unsteady lift forces. In this section, we are concerned with the inter- 
action of a blade row with the potential flow field of another blade row. 
The analysis is based on the formulas developed by Kemp and Sears 
(ref. 4). 

The following assumptions were made to compute the induced veloc- 
ities at an airfoil resulting from the potential flow field of all the 
blades of an adjacent blade row. 

• The flow through a rotor-stator combination is represented by 
nonviscous, incompressible flow through a two-dimensional cascade. 

• The airfoils of both blade rows, the velocity inducing and the 
lift producing, are represented by vortex sheets. 



The steady vorticity representing a blade is located on a straight 
line, that is parallel to the mean flow velocity through the 
cascade and that has the length of the local blade chord. 

Only the unsteady lift forces resulting from the relative motion 
of an airfoil through the potential flow field of the steady vort- 
icity of a blade row are considered. All unsteady lift forces 
resulting from secondary interactions are neglected. 

Only the transverse component of the induced velocity is considered 
in the computation of the unsteady lift. 


Four different interactions are possible between the potential 
flow fields of a rotor-stator combination: 

1) Rotor downstream of stator, unsteady lift induced at stator 

2) Rotor downstream of stator, unsteady lift induced at rotor 

3) Stator downstream of rotor, unsteady lift induced at rotor 

4) Stator downstream of rotor, unsteady lift induced at stator 

The analysis will be summarized for the first case, then the 

equations will be generalized so that they can be used for all four 
interactions. In reference 4 the following equation was derived for 
the induced, unsteady velocities at a stator blade located upstream of 
a moving rotor blade row: 
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with 
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temporal frequency 

complex frequency of the chordwise velocity 

distribution 

stagger angle 

chord length 

midchord plane separation 
mean relative flow velocity 
Glauert coefficients of n^ order 


40 



Index 1 = stator parameter 
Index 2 = rotor parameter 

The Glauert coefficients are the coefficients used in the series repre- 
sentation of the steady vorticity of an airfoil introduced by Glauert 
(ref. 38). 

The geometrical relationships upon which these equations are based 
are shown in figure 7. The unsteady, induced velocities described by 
equation (3.2.1) are characterized by two properties: 

1) They decay exponentially with increasing distance from the source. 

2) They are not medium bound and can, therefore, move at speeds 
different from the medium velocity. 

For these reasons, they cannot be treated as frozen convected gusts. 

Equation (3.2.1) defines the induced velocities only along the chord 

of the zero th stator blade. A small modification of the term g . will 

0,1 

make it possible to compute W(z,t) at any point of the stator flow field. 
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can be modified using the following approximation 


-tan = tan 


Y 2 - 


pM t 

%,1 COS Y 1 


(3.2.8) 





C 1 

g o, 1 “ 2p N 2 



H 


o,2 



M*. 


2 M 


M,2/J 


( 3 . 2 . 9 ) 


For the computation of the unsteady section lift we need to know 
the induced velocities at the midchord point. For the zero*'* 1 stator it 
is 
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The midchord plane is defined by the relationship 
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and the induced velocities at its location are defined as follows: 
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With the following definition of the blade separation distance, 
the induced velocity at the midchord of the .j blade can be computed 

I 

Y. 9 

J = _ 2 7T p . 

cos 

0.2.13) 


“j (t) 


TtC, 


E 

0=1 


... 2ir . . 
-xaN 2 jr J- 

Vl e 


lin 




( ) . 2 . 14 ) 

In reference 4 it is shown that for the transverse gust velocities 
described by equation (3.2.1), the lift response function assumes 
the following form: 
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where : 


J o ,J i 


H (2) h (2) 
H o ,H 1 


reduced frequency 

reduced complex frequency of chordwise velocity 
distribution 

Bessel functions of the first kind 
Hankel functions of the second kind 


Knowing the lift response function and the transverse gust velocity 
at the midchord point, the unsteady section lift can now be computed: 
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( 3 . 2 . 18 ) 


The Kemp-Sears lift response function is referred to as K^. If the 
imaginary part of X^ is equal to zero, and if the real part is equal 
to v^, then becomes the Sears function. If the imaginary part of 
X^ is not equal to zero, it increases proportionally to a and K^, 
then diverges for large values a because of the character of the 
Bessel functions. But, even then will the unsteady section lift assume 
finite values. This is due to the exponential decrease of the term 
g Q ^ with increasing values of a. It can be shown that the section lift 
remains finite for large values of o, if the following restriction is 
satisfied: 


\ (t^cos Y;l + c 2 cos y 2 J < d 


( 3 . 2 . 19 ) 


which indicates that the two interacting blade rows do not overlap, thus 
allowing relative motion. 



If the amount of turning is small in each of two interacting blade 
rows — that is, if the relative inflow and exit flow angles of a blade 
row can be set equal to the stagger angle — then can be approximated 
by the following equation: 
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(3.2.20) 


In this case, A-^ is, therefore, only a function of for a given 
combination of stagger angles. In figure 8, Kemp-Sears lift response 
functions are shown for three different combinations of stagger angles. 

Due to the convention (see appendix A) used for the Fourier transforma- 
tions, the complex conjugate of the section lift computed in equation 
(3.2.18) will be used in this analysis and the summation over 0 is 
modified to extend from -» to +°°. The new formula for the unsteady 
section lift is therefore 
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Inlet stator-rotor interaction 
unsteady lift response/ of rotor. 

= Rotor stagger angle 
T'lS = Inlet stator stagger angle - 45° 

7|s = Inlet stator stagger angle — 45° 

7r = 90° gives Sears function 

u = Reduced frequency 

FIGURE 8. - KEMP-SEARS LIFT RESPONSE FUNCTION 
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Since this is a lifting-line model, equations (3.0.2) and (2.2.17) 
will be used to compute the mode amplitudes. Therefore, the temporal 
Fourier transform of L^(p,t) has to be determined. 
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The duct mode amplitudes can now be computed using equation (2.2.17). 
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The equations for the other three potential flow field interactions 
possible in a rotor-stator combination are summarized in appendix F. 
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The computer subprogram AABAA (see section 3.1.2, volume II) computes the 
factor multiplying the frequency spike in equation (3.2.29) with to in the 
propagating region. 


3.3 Rotor-Alone Models 

A ducted fan will seldom, if ever, consist of a rotor alone, but the 
velocity disturbances at the rotor can be separated, in principle, into 
the component interactions and non-component-induced inflow distortions. 

The latter are the disturbances considered in these models. An exit stator 
will straighten out the flow so that upstream and downstream of the stage 
the flow will be axial. The inflow disturbances can be categorized as 
being either steady or nonsteady distortions. Steady distortion can 
occur under different conditions including for the fan stage of aircraft 
engines all conditions in which the air enters the duct inlet at a rel- 
ative angle of attack, and the distortions which result from the inlet 
and duct contours. Under static engine test conditions such distortions 
can result from ingesting the persistent nonuniformities in the air 
characteristic of the facility. Nonsteady distortions can be further 
classified according to their duration at the rotor. The classification 
parameter might be, for instance, the duration of the distortion times 
the blade passing frequency of the rotor. Steady distortions then lie 
at one end of this parametric scale, while small-scale turbulence lies 
at the other end. To have well-defined tone noise resulting from the 
distortion, the parameter must be, from purely intuitive reasoning, 
greater than unity. One possibility for the occurrence of distortions 
in this range is atmospheric turbulence, whose eddies can be stretched 
as they accelerate into the duct (for low-speed flight or ground static 
testing); see Pfenninger (ref. 18) and Hanson (ref. 19). It is a 
particular representation of this possibility that is developed into the 
model discussed in section 3.3.2. 

Figure 9 illustrates the rotor blade section velocity triangles. As in 
the component interaction models, the velocity disturbances must be 
determined in the coordinates attached to the blades. This is done in 




appendix G for an unspecified distortion, use of which is then made in 
sections 3.31 and 3.32 for the particular models. 

3.3.1 Steady Distortion . — When the distortion is expressed as the 
ratio of the deviation from the mean to the mean velocity, 
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Substituting this result into equation (2.3.1) gives, from equations 
(3.0.1) and (2.2.19), 
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The sum on j can be performed to give 
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The expression for the mode amplitude then becomes 
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and is given by equation (3.1.39) for this v^, is used to approx- 

imate V, i.e., y is taken as the mean stagger angle. All references to 
subscripts 1 and 2 have been dropped since there is only the one compo- 
nent under consideration, the rotor. The computer subprogram BCDAA 



computes the factor multiplying a frequency spike for w in the 
propagating region. 

This model is analytically similar in all respects but one to the 
inlet stator-rotor mode of the viscous wake interaction model. The 
exception is in the origin and specification of the velocity disturbances. 
In the viscous wake model, studies had been made, data had been correlated, 
and some theoretical ideas had been employed (Prandtl's mixing length 
theory) to determine the empirical formulas for the wake defects of air- 
foils (Silverstein, et al., ref. 37). Although much data has been gathered 
on the steady distortion of fan duct inlets, no work known to the 
authors has been published in which this data has been studied with the 
view of deriving a general empirical formula for computing the velocity 
distortion. In lieu of a sound theoretically or experimentally supported 
procedure, acousticians have begun to develop rather crude empirical 
models of their own in order to perform calculations. Wright (ref. 15), 
Lowson (ref. 17), Lowson and Ollerhead (ref. 39), and Barry and Moore 
(ref. 40) have investigated a "power law" for "blade loading harmonics." 
That is, in the present language, they represented the velocity distortion 
at a rotor disc in a Fourier series on the polar angle about the axis, 
at some average radius, and assumed that the coefficients raonotonically 
decreased as l with l the Fourier series index and q some positive 
number [loading harmonics are actually velocity harmonics times lift 
responses; see Hanson, chapter 6 of (ref. 19) for a discussion of this 
"law"]. This model is included in the computational possibilities of 
the computer subprogram BCDAA. It requires knowing besides the power, 
q, the first harmonic coefficient of the cosine series, the determination 
of which can possibly be determined from on axis microphone measurements 
(see Barry and Moore), i.e.. 



The other options in this subprogram include inputting the Fourier 
cosine and sine series coefficients as determined by the user (a simple 
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integration scheme is supplied outside the subprogram for computing these 
coefficients from distortion data), and a simple analytical model, the 
"cone model." This model was part of a study intended to clarify the 
power law discussed above and is given in terms of the formula 
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where : 

= ratio of maximum to minimum distortion 
A = span position of the maximum 


When the coefficients of the cosine and sine series are input, then 
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(3.3.10) 


and 



(3.3.11) 


with a^, the cosine and sine coefficients, respectively. 

3.3.2 Nonsteady Distortion . — Reference is made to Pfenninger (ref. 
18) for a discussion of the origin and occurrence of stretched eddies in 
fan ducts. Hulse, et al. (ref. 9) put forth the idea of isolated 
"patches" of distortion convecting through a rotor to explain some of 
their results, and Hanson (ref. 19) has gathered some data relevant to 
the occurrence of stretched eddies in fan ducts. This model puts into 
an analytical form some of these ideas. In appendix G an expression for 



the velocity perturbation on a rotor blade due to the convection through 
the rotor of a stretched eddy is developed; see equation (G23) . 


Substituting for u from this equation into equation (2.3.1), there 
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results for the dipole surface density function the expression 


2 ,=-“ 


SIN y W 




+ 


U) 


+ £M T 


(3.3.12) 


z£ 2 I M 


4,1 M 

+ cos T V 2 7 TT 


th 

where W is the £ coefficient in the complex Fourier series represen- 
tation of the axial component of the eddy cross section at the radius p g . 
The cross-sectional profile of an eddy is time independent, the assump- 
tion being that shape changes with z or z - M t are less significant 

^ th 

than magnitude changes. The term W is the £ coefficient in the ser- 
ies for the circumferential component of the eddy cross section at the 
radius p . Assuming the stretched eddy has cylindrical symmetry about 
its own axis, then a reasonable choice for a cross-sectional profile is 
the Gaussian (see fig. 10), 
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, ^ 2a z 
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z ^ s, sj z 

and similarly for W, (p ,d> ) with W, and a,, where W ,W, and a ,a, are 
the maxima of the distortions and the "transverse length scales” of the 
distortions, respectively. The location of the cross-sectional center 
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location of eddy center C at t = 0 

axial position of C referenced to the rotor midchord plane 

eddy length scales 

axial eddy velocity distribution 


FIGURE 10. - EDDY VELOCITY SPATIAL DISTRIBUTION AT TEMPORAL ORIGIN ( t = 0 ) 
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of the eddy is given by j R, $ J . This gives for w z ^(Pg) 
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with I the modified Bessel function of the first kind of order Z. 
Similarly, 
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Expressions for A and A are developed in appendix G, equation (G22), 
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where: 



T <p = L / M 
4 > 1 z 


= axial length scales of the eddy 

t = time at which the axial center of the eddy coincides 

with the midchord plane of the rotor 

Upon substituting these results into equation (3.0.1), the sum on 
j gives 
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with 1 = m-aN. Substituting this result into equation (2.2.19) for the 
mode amplitude gives 
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This spectrum for the mode amplitude is two series of Gaussians 
centered on the harmonics of blade passing frequency, with two-thirds of 
the power of each bump within plus or minus 2 tt/T z and 2ir/T^. For purposes 
of computing a single mode amplitude for each harmonic, the Gaussian 
"bumps" are replaced by weighted delta functions (see fig. ll).: 
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where 2B is the bandwidth of the filter used to define the tone and the 

subscripts z and <J> refer to the z and <f> components, respectively. The 

frequency-dependent parts of A" (oj) are then evaluated at harmonics of 
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with 


£ = m - a N 


Two methods are used to evaluate the chordwise integration of 
equation (3.3.24). One is to employ the pressure difference function 
of Naumann and Yeh (ref. 35) as was done in the steady distortion model 
of section 3.3.1, giving 
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with given by equation (3.1.39). The other method is to assume a 

lifting-line model with the lift response function developed by Filotas 
(ref. 36) for oblique gusts, the assumption being that the radial varia- 
tions of a compact eddy present to the blade oblique gusts made up of 
harmonic gusts in the chordwise and radial directions. Following in 
spirit the calculation "Distortion of Standing Wave by Strip" in 
chapter 11 of Morse and Feshbach (ref. 29), for calculating the section 
lift response of a three-dimensional gust, the blade is assumed to be a 
two-dimensional airfoil. The radial variation of the incident gust is 
decomposed into harmonics that travel up and down the span. These 
harmonics combined with the chordwise harmonic gusts produce oblique 
gusts at a given span position as required in the formulation of Filotas 
(see fig. 12). 
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FIGURE 12. - SCHEMATIC OF OBLIQUE HARMONIC GUST 


For this purpose the most significant radial variation in W „ and W,„ 

z x. 9 Jo 

is isolated, then the rectangular coordinate x is substituted for 
P g in this factor to derive the plane wave expansion. Since the most 
significant radial variation resides in the Gaussian factor of equations 
(3.3.15) and (3.3.16), the exporential and the modified Bessel function 
tend to be inversely related; then 
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These results produce oblique gusts of the form 
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(3.3.28) 
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with 


x-Vk 2 ^*) 2 


where, from appendix G, 


v 


£ 


c _£_ 
2 p s 


SIN Y 


and 



( 3 . 3 . 29 ) 


( 3 . 3 . 30 ) 


Since the gust wave form used by Filotas is with ^ -»■ -ip and since 
(see appendix H) 


T(h, - ip) = T + (h, i»). 


( 3 . 3 . 31 ) 


it is the complex conjugate of the oblique gust response function that 
is required here. Computing the section lift resulting from each k and 
adding them up then gives 
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Similarly, 


for 



(P g ,C",a)-nM T ^ d£' 


with W, and a, in place of W a . 
<f> 4> r z, z 


Substituting these results into equation (3.3.24) then gives the 
mode amplitude by this method. Subprogram BBCAA computes the mode 
amplitudes for w in the propagating region using one of two methods, 
i.e., the factor multiplying the frequency spike in equation (3.3.25) 
or that in equation (3.3.24) with the substitutions from equation 
(3.3.32) and the w equivalent equation. 


3.4 Numerical Considerations 

Each of the models produce an expression for a mode amplitude which 
can be written quite generally in the form 

1 


n 

(3.4.1) 

where a is the harmonic index specifying the frequency under considera- 
tion. To compute all of the upstream- or downstream-propagating mode 
amplitudes for this harmonic, the set of m's and n's for these modes 
must be determined from the cut-on criterion. 
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where : 

N = number of rotor blades 
K 

E = a bound for the eigenvalues to be computed . 

.D 

Having determined the indexes for the propagating modes, the corresponding 
set of y * s are computed. These preliminaries determine the set of 
integrals that have to be computed. The integrand of each of the inte- 
grals is then subjected to the following factorization: the factors that 

are independent of the radial, or span, variable p , i.e., the constant 
factors, are identified and placed outside the integral sign. Then the 
integration interval is subdivided into a number of equal subintervals, 
the number of subdivisions being determined by the number of oscillations 
of the most oscillatory factor. The remaining p s ~dependent integrand is 
then factored into a part that varies significantly, possibly changing 
sign, within the subinterval, and a part that varies slowly, allowing its 
average value to represent it within the subinterval. The calculation 
can then be written symbolically in the form 


A± 

mno 


N 


= (CONSTANT) 



j=l 


OF SLOWLY VARYING FACTOR) 


3 j+l 


(OSCILLATORY FACTOR) dp s 


( 3 . 4 . 3 ) 


with a^. the left endpoint of the ^ th subinterval. The integration on 
each subinterval is then performed using an 8-point Gaussian rule. 
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U.o USERS GUIDE TO COMPUTER SUBPROGRAMS 


These subprograms compute all the mode amplitudes for a given harmon- 
ic of blade passing frequency for either upstream- or downstream-propagating 
modes. Subroutine AAAAA does this for the viscous wake interaction model, 
subroutine AABAA for the potential field interaction model, subroutine 
BCDAA for the rotor in steady distortion model, and subroutine BBCAA for 
the rotor in nonsteady distortion model (see fig. 13). 


4.1 Calling Sequence 


All four primary subroutines have the same calling sequence and 
are called as follows : 


DIMENSION MUSE (MDIM) ,MAXN (MDIM) , 

* ARMUMN(NDIM,MDIM) ,ARMISC(40) , AR (MAXDIM , MAX J , 3 ) 

COMPLEX ALP HAMN(NDIM, MDIM) 


CALL < 


AAAAA 

AABAA 

BCDAA 

BBCAA 


* 


( ARMI S C , MAXDIM , MAXJ , AR , MDIM , NDIM , 


ARMUMN , NOFM , MUSE , ALPHAMN , TERROR) 


4.2 Input, Input — Output, Output 


INPUT 

ARMISC an array of dimension 40 used for input and fully 

described in the section below. 



AAAAA — Viscous wake 
interaction 

AABAA — Potential field 
interaction 

BCDAA — Rotor in steady 
distortion 

BBCAA — Rotor in 
non-steady 
distortion 


A B 



FIGURE 13. - PRIMARY SUBROUTINES 












MAXDIM 


HAXJ 


The first variable dimension of the array AR, which 
corresponds to the radial positions and must be 
dimensioned at least 2 more than the largest number 
of radial positions to be input for any component 
(inlet stator, rotor, or outlet stator). 

The second variable dimension of the array AR, cor- 
responding to the number of aerodynamic and geometric 
variables input and, per package, must be at least 


Package 

1 (AAAAA) 

2 (AABAA) 


3 (BCDAA) 


4 (BBCAA) 


Value 


11 


9 + max | aRMISC (I) : 1=18,19, 20 J 

where : 

max ARMISC(I) _> 3, 1=18,19,20 

11 if ARMI SC ( 2 2 ) = 0,1 

12 if ARMISC(22) = 2 
11+ARMISC(23) if ARMISC(22) = 3 

9 if ARMISC(25) = 4 
11 if ARMISC (25) = 3 


AR 


A three-dimensional array of geometric and aerodynamic 
variables, described in the section 4.4. 


MDIM 


NDIM 


The variable dimension corresponding to the maximum 
number of spinning mode indexes. MDIM should be at 
least 2Eg + 1 but need be no more than 201 to ensure 
that all spinning mode indexes can be calculated 
(E = w*; see equation [2.1.8]). 

D 


The variable dimension corresponding to the number of 
radial mode indexes, limited to at most 40. 



INPUT/OUTPUT 


ARMUMN 

NOFM 

MUSE 

MAXN 

OUTPUT 

ALPHAMN 


The real array-dimensioned NDIM x MDIM-that contains 
the eigenvalues calculated in the primary subroutine. 
Upon subsequent calls to the primary subroutine the 
eigenvalues from the previous case are reused if the 
parameters used in calculating the eigenvalues remain 
unchanged. 

The number of spinning mode indexes computed by the 
program; can be reused as discussed in ARMUMN. 

An integer array-dimensioned MDIM-that contains the 
NOFM spinning mode indexes and can be reused as 
discussed under ARMUMN. 

An integer array-dimensioned MDIM- containing the 
maximum radial mode indexes corresponding to the 
spinning mode indexes in array MUSE. Array MAXN can 
be reused as discussed in ARMUMN. 


The complex array-dimensioned NDIM x MDIM-containing 
the calculated mode amplitudes. For any index IOFM 
( IOFM = 1..., NOFM), M = MUSE (IOFM) is an available 
spinning mode index. For this M, the available radial 
mode indexes are N = 1, 2,..., MAXN (IOFM). Then the 
corresponding mode amplitude is ALPHAMN (N, IOFM) and 
the corresponding eigenvalue is ARMUMN (N, IOFM) . 



ERROR RETURN 


IERROR 0 = all calculations have been complete successfully 

2 = an incomplete set of eigenvalues is calculated, 

and the corresponding mode amplitudes are calculated 

4 = there are no eigenvalues, hence, there are no mode 
amplitudes 


4.3 Input Array ARMISC 

In this and following sections, the code given below is used to refer 
to the subroutine packages : 


Package Code 

AAAAA 1 
AABAA 2 
BCDAA 3 
BBCAA 4 


In the following definitions, all variables referred to as non- 
dimensional are nondimensionalized as prescribed in appendix A. 


ARMISC 

ARMISC 

array Packages 

index used in Definition 


Nondimens ional average distance be- 
tween the midchord planes of the 
inlet-stator and the rotor; see 
figure 4. 


1 


1 , 2 



ARMISC 

array Packages 

index used in Definition 


2 1,2 


3 1,2, 3, 4 

4 1 , 2 , 3, 4 

5 1,2 


6 1 , 2 , 3, 4 


7 1 , 2 , 3, 4 

8 1,2 

9 1,2 


Nondimens ional average distance 
between the midchord planes of the 
rotor and the outlet-stator vanes; 
see figure 5. 

The hub-to-tip ratio of the duct 

-1 for upstream sound propagation 
1 for downstream sound propagation 

1 indicates inlet stator-rotor inter- 
action 

2 indicates rotor-outlet stator inter- 
action 

Trace printout option, where 

0 = no trace printout 

1 = trace printout of major factors 

in the primary subroutine 

2 = detailed printout of the eigen- 

value calculation in subroutine 
ZEROS; see section 3.2.2 of 
volume II 

Mp the rotor blade tip Mach number 
^ISV’ t ^ e number of inlet-stator vanes 
Nqsv> the number of outlet-stator 


vanes 



ARMISC 

array 

Packages 


index 

us ed in 

Definition 

10 

1,2, 3, 4 

N , the number of rotor blades 
Kd 

11 

— 

Not used 

12 

1 

4>Qg, in radians; phase angle for ad- 
justment of skewness of the incident 
wakes on the outlet stator; see 
appendix E. 

13 

1 

<j> , in radians; phase angle for ad- 
justment of skewness of the incident 
wakes on the rotor; see appendix E. 

14 

1,2, 3,4 

The harmonic index 

15 

1,2 

Z^g, the axial position of the inlet 
stator 

16 

1,2 

Zq£, the axial position of the outle; 
stator 

17 

1,2, 3, 4 

Z , the axial position of the rotor 
K 

18 

2 

-1 means the upstream blade row is 
the sound generator in a potential 
flow field interaction 
1 means the downstream blade row is 
the sound generator in a potentia 
flow field interaction 



Definition 


ARMISC 

array Packages 

index used in 


19 2 


20 2 


21 2 


22 3 


The number of inlet stator vane 
Glauert coefficients, which is 
0 if ARMISC (5) = 2 or 
ARMISC (18) = -1 

n where 3 <_ n 5. 15 
if ARMISC (5) = 1 and 
ARMISC(18) = 1 
(see sec. 3.2) 

The number of rotor blade Glauert 
coefficients, which is 

0 if ARMISC (5) + ARMISC (18) - 1 or 2 

n where 3 5 n < 15 

if ARMISC (5) + ARMISC(18) = 0 or 3 
(see sec. 3.2) 

The number of outlet stator Glauert 
coefficients, which is 
0 if ARMISC (5) = 1 or 
ARMISC (18) = 1 


n where 3 <_ n 15 
if ARMISC (5) = 2 and 
ARMISC (18) = -1 

(see sec. 3.2) 

0 = no distortion 

1 = distortion is represented by the 

cone model 
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ARMISC 


Packages 
used in 


Definition 


array 
index 

2 = distortion is represented by the 

power model 

3 = distortion coefficients are input 

23 3 V a/ V i if ARMISC ( 22 ) = 1 

q if ARMISC (22) = 2 

Total number of distortion coeffi- 
cients if ARMISC (22) = 3 

24 3 A if ARMISC (22) = 1 

Skip factor used with distortion coef- 
ficients if ARMISC (23) = 3 

25 1,2, 3, 4 Defines the lift function to be used: 

2 = LIFTFN2, the generalized Sears 

lift response function; see equa- 
tion (24) of appendix I; used in 
package 2 only. 

3 = LIFTFN3 or NONCPT; LIFTFN3 is the 

combination of lift response func- 
tions as developed in reference 6; 
NONCPT is the lift response func- 
tion for noncompact source theory, 
see equation (22) of appendix I; 
used in packages 1, 3, and 4. 

4 = LIFTFN4, Filotas lift response 

function; see equation (25) of 
appendix I; used in package 4. 
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ARMISC 

array Package 

index used in Definition 


26 

27 

28 4 

29 4 


Not used 
Not used 

R, radial position of the eddy center 

4>, angular position of the eddy center 
in radians 


30 4 


31 4 


32 4 


33 4 


W z> axial eddy velocity component at 
the eddy center, nondimensionalized 
with the average axial flow velocity; 
see figure 10. 

W^, angular eddy velocity component 
at the eddy center, nondimensionalized 
with the average axial flow velocity; 
see figure 10. 

a 2 » eddy length scale in the direc- 
tion normal to the average flow 
velocity for the axial eddy velocity 
component; see figure 10. 

a , eddy length scale in the direc- 
tion normal to the average flow 
velocity for the angular eddy velocity 
component; see figure 10. 



ARMISC 

array 

index 

34 

35 

36 

37 

38 

Note: 

39 

40 


Packages 

used in Definition 


4 L^» eddy length scale in the direc- 

tion of the average flow velocity 
for the axial eddy velocity compo- 
nent; see figure 10. 

4 L , eddy length scale in the direc- 

<P 

tion of the average flow velocity for 
the angular eddy velocity component; 
see figure 10. 

4 B, upper bound of the frequency band 

considered in the generation of' tone 
duct mode amplitudes by nonsteady 
distortion; see figure 11. 

4 t, time when eddy center is located 

in rotor plane 

1,3,4 0 compact source (LIFTFN3 is used) 

#0 noncompact source (NONCPT is used) 

ARMISC(38) can be used if and only if ARMISC(25) = 3. 

Not used 

Not used 


4 ,4 Input Array AR 


AR is an array of dimension MAXDIM x MAXJ x 3 which contains geometric 
and aerodynamic data as either average values across the span or as values 



given at spanwise positions. The values of the array AR, AR(I,J,K), are 
defined according to the use of each index I, J, and K. 


AR 



array 

Package 


index 

used in 

Definition 

T — ( 

II 

1,2 

Refers to inlet stator data and is 



used when ARMISC(5) = 1. 

K = 2 

1, 2,3,4 

Refers to rotor data, i.e.. 



AR(l,J,2) contains the rotor data. 

K = 3 

1,2 

Refers to outlet stator data and is 


used when ARMISC(5) = 2. 

0 if average value 
of quantity cor- 
responding to J 
and K indexes is 
to be used 

n(K) if spanwise data 
corresponding to 
J and K indexes 
is to be used 
where n(K) = num- 
ber of spanwise 
positions for 
index K 

1=2 1,2, 3,4 AR(2,J,K) refers to the average value 

of the quantity corresponding to J 
and K indexes . 
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AR 



array 

Package 


index 

used in 

Defini tion 

I = 3,4. .. 

1,2, 3, 4 

AR(I,J,K), I = 3,4 ... refers to 


spanwise data corresponding to J and 
K indexes. AR(3,J,K) refers to first 
nondimensional duct radial position, 
AR(4,J,K) refers to second nondimen- 
sional duct radial position, and so 
on, in monotonic increasing order. 

1,2, 3, 4 p, nondimensional duct radial position 

1,2, 3,4 C, nondimensional chord 

Not used 

1 C^, drag coefficient 

Not used 

1,2,3, 4 dC^yda, derivative of steady-state 

lift coefficient, C^, with respect to 

incident angle, a 

J = 7 1,2, 3, 4 M^, relative inflow Mach number of a 

blade row; see figure 3. 

J = 8 1,2, 3,4 M , relative exit flow Mach number of 

£ 

a blade row; see figure 3. 

J = 9 1,2, 3, 4 M , axial flow Mach number; see fig- 

ure 3. Note: the average value, 

AR(2,9,K), must always be given. 


J = 1 
J = 2 
J = 3 
J = 4 
J = 5 
J = 6 
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For the remaining definitions of J, the subroutine packages are dis 
cussed separately. 

For package 1 (AAAAA) : 

AR 

array 

index Definition 

J = 10 f, the ratio of maximum blade camber to the half 

chord 

J = 11 a, the blade angle of attack 

For package 2 (AABAA) 

Let NGC = max JaRMISC(18+K) : K = 1,2,3} 

= number of Glauert coefficients (see sec. 3.2) 


AR 

array 

index Definition 

J = 10 A°, Glauert coefficient of order 0 


J = 11 



Glauert coefficient of order 1 


J = 9 + NGC 


NGC-1 

A 


Glauert coefficient of order NGC-1 



For package 3 (BCDAA) : 


If ARMISC(22) =0, 1: 

AR 

array 

index Def ini tion 

J = 10 f, the ratio of maximum blade camber to the 

half-chord 

J = 11 a, the blade angle of attack 

If ARMISC(22) = 2: 

AR 

array 

index Definition 

J = 10 f, the ratio of maximum blade camber to the 

half-chord 

J = 11 a, the blade angle of attack 

J = 12 which is used in the power model; see 

equation (3.3.8) 

If ARMISCC22) = 3: 

Let NDC = ARMISC(23) = total number of distortion coefficients and 
let SF = ARMISC(24) = skip factor. 
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Definition 


AR 

array 

index 


J = 10 


f, the ratio of maximum blade camber to the 
ha If -chord 


J = 11 

a, the blade angle of attack 

J = 12 

a^p, distortion coefficient 

J = 13 

bgp, distortion coefficient 

J = 14 

a 2*SF’ distortion coefficient 

J = 15 

• 

^2*SF* distortion coefficient 
• 

• 

• 

J = 10 + NDC 

• 

a NDC » distortion coefficient 

2 * SF 

J - 11 + NDC 

b , distortion coefficient 

2 * SF 

For package 4 (BBCAA) : 


If ARMISC(25) = 4, 

no further J's are required and ; 

complete; if ARMISC(25) 

= 3: 



AR 

array- 

index Definition 


J = 10 f, the ratio of maximum blade camber to the 

half-chord 


J - 11 


a, the blade angle of attack 


4.5 Restrictions and Limitations 

The use and restrictions on the input arrays ARMISC and AR and the 
special input/output NOFM, MUSE, MAXN, and ARMUMN are given in the 
previous section. 

The maximum spinning mode index is limited (see subroutine EGHVAL2, 
sec. 3.2.1 of vol. II) in absolute value to 100 and the maximum radial 
mode index is at most 34. 


4 . 6 Diagnos tics 

Diagnostic messages related to the computation of eigenvalues are 
provided according to the printout control parameter ARMISC(6) and the 
error return parameter IERROR. 

When ARMISC (6) j4 0 and IERROR = 2, the following is printed: 

A REDUCED SET OF EIGENVALUES IS AVAILABLE 
COMPUTATIONS WILL PROCEED. 

and when ARMISC (6) 4 0 and IERROR = 4, the following is printed: 

THERE ARE NO PROPAGATING RADIAL MODES 
NO COMPUTATIONS CAN BE MADE. 



4.7 Storage 


Each subroutine package with the subprograms listed in section 4.9 
requires the following approximate storage (octal) : 

Package Storage 


1 

(AAAAA) 

12,000 

2 

(AABAA) 

15,000 

3 

(BCDAA) 

12,100 

4 

(BBCAA) 

20,100 


4.8 Timing 

The timing in general is dominated by the calculation of the eigen- 
values. For the sample cases run, the average time per case is: 

Subroutine Time in decimal 

Package seconds 

1 55 

2 62 

3 117 

4 145 
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U.9 Internal Subprogram Relationships 

The following listing gives the subprograms required for each package 
except the standard CDC - issued system routines such as SORT, SIN, etc. 
An "X" indicates that the subprogram is used and a blank indicates 
nonuse. The subprograms used are described either in this documentation 
(vol. II), or NASA Langley Research Center library subprograms (refs. 

41, 42, 43), which are marked by an asterisk. 


Subprogram 


AAAAA AABAA BCDAA BBCAA 


EGNVAL2 

ZEROS 

EQATION 

UNEGNFN 

EGNNORM 

FACTINT 

FACTIN2 

FACTIN3 

FACTIN4 

LIFTFN2 

LIFTFN3 

LIFTFN4 

DISINT 

FUNIN4 

NONCPT 

APR0X1 

APR0X2 

JARRATT 

GAUSS 

GAUSS 2 

BSSLS 

BESNX 

BESJLA 

BESIE 

BESIK 


X 

X 

X 

X 

X 

X 


X 

X 

X 

X 

X 

X 


X 

X 

X 

X 

X 


X 


X 


X 


X 


X 

X 

X 

X 

X 

X 

X 


X 

X 

X 

X 

X 


X 

X 

X 

X 

X 

X 

X 

X 

X 


X 

X 

X 

X 

X 


X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 



Subprogram 

AAAAA 

AABAA 

BCDAA 

BBCAA 

ROCABES 


X 


X 

SICI 




X 

GRTHFCN 




X 

BF4F* 

X 

X 

X 

X 

MTLUP* 

X 

X 

X 

X 

ALGAMF* 


X 


X 



5.0 RESULTS AND CONCLUSIONS 


Analytic models for computing the sound pressure at harmonics of the 
blade passing frequency generated by a single stage of a fan or 
compressor operating in an infinite, hardwall annular duct have been 
developed and programmed to give numerical results on a CDC 6600 
computer. The outputs of the programs are the hardwall, plug flow mode 
amplitudes for those harmonics above cutoff— the propagating-mode 
amplitudes. Using the mode function of section 2.2, the pressures at 
field points in the duct on either side of the fan stage can be com- 
puted by first adding up the product of the mode amplitudes times mode 
functions to get p and then using equation (2.2.4) to get the pressure. 
Since the pressure is nondimens ional, the appropriate reference pressure 
will have to be used to get the correct SPL level. A simple procedure 
is to compute the SPL from the computed pressure and add 197 dB, the SPL 
level of standard atmospheric pressure. However, the more important use 
to which the results of the programs can be put is as inputs to a duct 
acoustic program such as the one envisioned by Zorumski (ref. 24). 

No attempt has been made to seek out data to compare with the results 
of this analysis. This should be done, but the data should be well 
gathered. That is, the experiments should be performed under sufficiently 
controlled conditions so that a determination can be made on which part 
of the models— the incident velocity disturbances or the acoustic response 
of the blade rows— gives the variance. 

The interaction models produce the rotor-stator blade number ratio results 
of Tyler and Sofrin, and Lowson, and they allow studies to be made of the 
relative importance for sound generation of the potential fields versus 
the viscous wakes and the fall-off of sound level with blade row spacing. 
Among other possibilities, the rotor-alone models could be used to gain 
a better understanding of the differences between static test data and 
flight data by considering the different inflow conditions of the two 
situations . 



It is felt that to improve these models, one would seek better definition 
of the velocity disturbances and a less restrictive approximation to the 
acoustic response function involving, perhaps, the carrying through of the 
formulation of appendix D. 


Boeing Commercial Airplane Company 
P.0. Box 3707 


Seattle, Washington 98124, May 31, 1974. 



APPENDIX A 


In this report, the temporal Fourier transform of a real function f(t) 
is referred to as the spectral density of f and is defined by 


co 

/■ 


ltilt 


f (to) = / f(t) e 


dt 


so that 


(Al) 


f(t) 


1_ 

2tt 



— icot 
e 


_ CO 


dco 


(A2) 


A function that is periodic has a spectral density that is a sum 
of Dirac delta functions, or frequency spikes, i.e., if 


f(t) 



f n e 


■ioi t 
n 


then 

00 

f ^ = f n 5 

n=-» 


(A3) 


(A4) 


All complex functions will obey a symmetry principle in m (when 
so that only positive frequencies need to be considered in the 
final calculations, i.e,. 



CO 


f(t) = J 2 Real ^f(w)e iwt ^ du> 
o 

Spatial Fourier transforms are defined by 


(A5) 


g(K) 


00 

J g(x> e~ iKx dx 

— oo 


(A6) 


with 


00 

g (x) f i( K > eiKXdK 

(A7) 


and Fourier series, such as on the interval 0 to 2ir, 


2ir 

/ -i£<J> 

Vl(<f>) e “(f* 

o (A8) 


with 


OO 

w( *> - b £ w * 

£=-oo 


(A9) 
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The linearized differential equations of an inviscid, compressible 
fluid in which isentropic perturbations from equilibrium can take place 
are assumed to be sufficient for the analysis. The speed of sound is 
then taken to be a constant, with 


speed of sound 



(A10) 


where p is the pressure perturbation and 6 is the density perturbation. 
It is assumed that for computing the perturbation variables the mean 
velocity can be taken uniform and constant. Since everything will take 
place in a duct, the duct outer radius will occur often as a factor 
multiplying an inverse length type quantity. 

The choice is made then to nondimensionalize all quantities from 
the beginning, choosing the duct outer radius, the speed of sound, and 
the mean fluid density as the basic scale factors. Letting these be, 
respectively, r Q , a Q , and d Q , and letting primed quantities be dimen- 
sional, then the nondimensional quantities of interest are 


6 = 



density perturbation 


(All) 


P = 



2 

a 

o 


pressure perturbation 


(A12) 


v = v' / 




- *7 
= ' r 


velocity perturbation (A13) 


mean velocity (A14) 


x 


length 


(A15) 



time 


(A16) 



o 


j* 

U3 = u)' o/a Q frequency (A17) 

The nondimensionalized basic equations are 


Dv 

Dt 


+ 



P = 


0 


(A18) 


— + "T * "v* = 0 

Dt 


(A19) 


6 = p 


with 


D_ 

Dt 


9 

9 t 


+ M 


(A20) 


(A21) 


and M is a constant vector. For computing sound pressures, M is taken to 
be axial and positive in the positive z direction, where the z-axis of 
the coordinate system is coincident with the duct axis. Substituting 
for 6 in equation (A19) and decoupling v from p by applying V* to equa- 
tion (A18) and D/Dt to equation (A19) , the convective wave equation for 
p is arrived at : 


-□ c p 



0 


(A22) 


This is the basic, nondimensionalized equation for this analysis. 



APPENDIX B 


The required particular integral to the inhomogenous equation, equation 
(2.1.10), is the familiar Helmholtz equation Green's function for out- 
going waves : 


iw*R* 
e 

4tt R* (Bl) 

with 


R* = "\j p 2 +p o 2 -2pp o COS ( (j)-^ ) + ( Z*-Z* ) 2 

(B2) 

It is more convenient in cylindrical coordinates to express this 
solution in a combined Fourier series-Fourier transform form. 


FF 





j »( w <K s ( z *" z °) ds 

(B3) 


(see, e.g., Levine and Schwinger, ref. 44) 


V 





and 


( p IF p> P Q 

( P D IF P< P 0 


(B4) 


(B5) 
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(B5) 


( p IF p< p 
P< = { 

Ip IF p> p 
l o o 

and H ^ is the Hankel function of the first kind, of order m, and 
m 

J is the Bessel function of order m. By separation of variables, the 
m 

general solution to the homogenous Helmholtz equation in cylindrical 
coordinates can be similarly expressed. 


CO 

ni=-a> _cc 


is(Z-Z*) 
e o ds 


(B6) 


with Y the Neumann function of order m, with A and B undetermined, 
m mm 

and the subscript on y n intended to convey the use to which the 

general solution will be put, i.e., as a boundary effect. It is by 

adding y „ to y ™ and solving for A and B from the boundary conditions 
B FF m ^m 

that the boundary effects are included in y. Adding them gives 


Y a 


Y + Y 
FF B 


” im ( ip-<p ) ( 

t 2 > ° J Vm (pp) + B m V PP > 

m=-» - co ' 


H ( ^(p P> ) J (up 
m * m ^ i 


* e 


la < Z *<> ds 


(B7) 


Employing the boundary conditions of equation (2.1.11), 


r — — = o AT p = l , P=n , 
3 P 
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gives the simultaneous equations for and B^, 


K J ’(y) 

m m 


,(D. 


+ B Y ' (y) « J (yp ) H v '(y) 
mm mom 


A 

m 


J'(yn) + B Y'(yn) = J'(yn) H^(yp ) 
m mm m mo 


(B8) 


where the prime denotes differentiation with respect to the argument. 
The solutions to these equations are 


A 

m 


= J (yp ) + i Y'(y) 
mo m 


Y '(yn) - Y (yp ) J'(yn) 
mom mom 

A (y) 
m 


(B9) 


and 


B = 
m 


-i J' (yn) 
m 


J m (w o ) Y m (|i) ' V’V J m (p) 


(BIO) 


with 


A (y) = J'(y) Y ' (yn) - J'(yn) Y' (y) 
m mm mm 


(Bll) 


Substituting these results into equation (B7) and performing the 
algebra gives for y , 


Y = 


2it 


im ( <p~<p j . / \ 

22 e v ° } M p ’ p °’ z *' z * ,a) *) 


m =_co 


(B12) 



where 


I 

m 


CO 

1 f Xm(y) is ( z -Z o) 

4 J zrury* 


(B13) 


and 


m 


■ < J (hP>) Y'(y) - Y (yp ) J' (y)i | J (yp ) Y'(yn)- Y (yp ) J'(yn)] 
1 m ? m m > mi i m < m m m i 


(B14) 


^ and A^, when analytically continued into the complex s-plane, are 
regular functions of the complex variable s. This can be seen by con- 
sidering A m and letting m* have a small positive imaginary part; 


+ i e y £ > o. 
K 


(B15) 


Then the branch points of y are raised off the real s-axis, as in fig- 
ure lU. The phases of y are indicated in this figure as well. Across 
the cut for y, the following relations hold; 


J m (y + ) = (-D m Jm(y") 


Y m (/) - C-l>" [ YJp-) + 2 i J m (lO _ 


CB16) 


Then, using the relationships 


2 J*(y) = J . (y) 
m m— 1 


^m+1 


(M) 


(317) 



(CJ* R + ie ) 


X 


REAL S 



FIGURE 14. - COMPLEX S PLANE 



(B17) 


2 Y’ (y) = Y (u) - Y (u) , 
m m-1 m+1 

it follows that 


A (y + ) = 4 
m 


(j 


m-l 


[^Vl^^l^) ( Vl ( - +) - Y m+l (w+)) 

(m+) - Vl (m+)) <Y m-l (,,+> - Vu (u+>> 1 


* > 
m 


(B18) 


from equation (B16) . Similarly for x^* 

Then, since the zeros of A are not coincident with the zeros of 

m 

the integrand of the integral for 1^ is meromorphic in the finite 

complex s-plane. For z* - z* > 0, the contour can be closed in the 

upper half-plane to encircle the poles at the zeros of A m * If these 

zeros (an infinite number, since A is transcendental) occur at s , 

m mn 

n=0, 1, 2, ...» then 


I 

m 


CO 




RES 


n=o 


( s-s ) 
mn 


X m is(Z*-Z*) 

A e 

m 


s->- s 

m 


(B19) 


The zeros in the lower half-plane occur at s = -s , n = 0, 1, 2, 

mn 

(see fig. 15), so that when z* - z* < 0, and the contour is closed 
in the lower half -plane. 



I = + — 
m 2 


(s + s ) !■ is(2*-zj)' 

v an ; a e 


s— -s (B20) 

ran 


Since A is a function of y, the zeros of A are y=y , n = 0, 1, 2, . .., 
ra m mn 

real, nonnegative numbers, and since 


A m = - A M = U 
m u m » mn -mn 


This gives for s , 
mn 


/ j.2 2 

s = \ w* - y 
mn V mn 


(B21) 


See appendix J for a discussion of the numerical scheme used to 
calculate these zeros. 


Near a zero in the upper half -plane. 


A (y) s~ A (y )+( s -s ) — A’ (y ) 

m m mn mn ds I m mn 

s = s 

ran 


■ < s - s > T- 

mn ds 


A 'm (M mn ) ’ 


(B22) 


x (U ) i s (Z*-Z*) 
m S5 e mn o 

ini a' ( y__) 

ds m mn 
s=s 

mn 


(B23) 
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•prc*ft 




Also, 


m 


( p \ = Y' ( p \y' (r\\i \ R ( p p^ R (p p\ 
y mn f m y mn y m I mn 1 m y mn J my inn °f 


(B28) 


so that 


RES 


Y ' fu Y' frill "V R ( p pN R (m p^~V 

m y mny m V. mny m V nm y m vmn Qy e 


is 


mn 


( Z *- Z o) 


^ 1 C 1 - 7 - 2 ) C"'»n) R » Ow) - n ( L - C>Vn)\ ow! 

mn ( M mn ^mn' 


!) 


(B29) 


Dividing top and bottom of this expression by Y^ (^ mn ^) ^ (j r '^ J mn ) 


gives 


RES = 


R (v P ) 
m V mn s 

R m C' J mn P o3 


is (z*-Z*) 

e mn u 

1 ^ (l * 

\ R m C^mn) / 

m 2 \ R m C nP mn) 

} smn 

y mn V p 2 

mn 

) y ; e»r \ 

nV ! V, OVD 

mn 

1 


(B30) 


Then, 


R fp ") 
mV mn2 


R 

m 


C P mn) 


t' fu h y' fu \ j fu ^-j' Cm \ y Cu V 

m V "in; mV. mny m V. mn^ m V. mny m V mn/ 


~ p R /p 'N 

2 mn ra V. ’mn ) 


(B31) 



and 


^m (^^mn) m tt 
Y m (nUrnn)” 2 Pmn 0 



m 



(B32) 


Thus , 


RES. = - 


(\i p') R (u p'N 

\ mn ' m V mn o/ 


xs 


m \ mn • m v. mn oj 

i { ( l - f) R - ^ • 1,2 1 1 ' ^r) R " <*4 

' u mn ' ' H mn 


mn 


& 


•Z*\ 
o ' 


mn 



<R f\i p's <R (\i p 
m v mn J m v - mn H o>' 


is ( Z*-Z*^ 
e mn v o J 


s 

mn 


(B33) 


where 


<R ( M p") 
m V mn J 



(B34) 


and 


N 

mn 


■ 7 { ( i) < ( ‘ ^ ( ^nfc) R “ K ) | 

mn 


(B35) 
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is the normaliz- 


This is a convenient form for the result since N 

mn 

ation of the R (y p) function, i.e., 
m mn 


/ «„(%,“) fl m(v) pdp ' S n l 
n 


Then, for Z* - Z* > 0, 


I 

m 


E d? / y <R /y p\ 

m l H mn H J m l H mn K o J 


e 


is (Z*-!*') 
mn V o / 

JTs 


mn 


For 


For Z* - Z* < 0, 

o 


dy 

ds 



s 

mn 



So that, following the same algebra, from equation (B20), 


I 

m 


E fl fp pU /li p] 

m l mn / ml mn o f 


-is (Z* - Z*) 
mn ^ os 


2is 


mn 


(B36) 


(B37) 


(B38) 


(B39) 



I'i n; i I I y , the result lor I is 

m 


I 

in 



7* - /.*, w*| = V'' A? p \ <U ( \i (. \ d ( 7* 

o f / a in l inn / I inn o / inn \ o / 

n=o 


( li A r i > 


wi th 


d mn •(Z*-Z*,u>*) = 


i s 

e mn 


7 .*~ 7 * 

o 


2is 


mn 


(BA1) 


Substituting this result into equation (1412) gives the final result, 
the solution for y: 



Z*,ID* i 
o J 

( 1142 ) 


Considering w* to be a complex variable, then s as a function of 

mn 

the complex w* has the branch line and phases as indicated in figure 15. 


Hence, for w* > g , s (-w*) = -s (w*), and when w* < p , s 

mn mn mn mn mn = 

V 2 2 

- w* when to* approaches the real axis from above it. The 

region of the real axis outside the branch points at ±p is referred to 

mn 

as the propagation (or "cut-on") region. It is seen that 



r\ mi 


mn 


(z* - Z*,a>*) 


(B43) 


where t means complex conjugation. This easures that y (-to*) = y 
Finally, just above the real axis in the propagating region, s has a 
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< >- 


positive imaginary part greater than that of w* so that for | z* - z*| -* 00 , 


„ -a Inn (u>*) |z* - Z*| 
e o 


(B4A) 
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APPENDIX C 


The evaluation of equation (2.2.7) is done using the transforms of F 

and D . Letting 
mn 



-icot 

o 

e dm 


and 


CO 



with t = t - t , then 
o 


(Cl) 


(C2) 


p(r^ t) = - 7T- S > S ' <R Ai e 1 " 1 ^ 

H 2tt / j / i mv®n J m v mn sy 

m=-°° n=o 


(C3) 


OO 

f f dm dm’ fm v ± 1 " , N irn’t 

*/ I — e , + K e F (w) e 

J J 2tt 2tt ^ p <f> mn z I 


2ir 


C z - z s-“D 2 / / 6 C*o - ♦ + e ° e 


-ira<j) -i(a>-u) , )t 


° dt 


£=— CD -OO O 



Considering the t and 4 > q integration, the delta function gives 


E 

I ” — rr 



- $ + 2-nSL-nt 




i(o)-w' )t Q -im<j> 0 
e e dt 0 d<}i o 


(C4) 


OO r,_j „ to— 01 v 

-Si- -) 

£=-oo 


Q 




'V .u-oj'+mn 


4> 

d(f 0 


o 



Q> -4> + 2ttjQ -im(j) 0 

e d4> 0 


Then, using Poisson's sum rule, and doing the <J> q integration, gives 


R.H.S . 



6 J)e^ CT ^ — 


in 


+ mn 


( 


-2iri 

e 


to-ui* + mn \ 

n _! J 


= (aQ-(u)-a)' j) e°^ 

0=— oo 


i 

a + m 


(• 


-2iri(a + m) 


-) 


= 6 (^an-(o3-aj' )} e ^2ir 6 a ,-m 

a=-® (C5) 
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271 5 Coj-tu* + nQ) e -im$ 


*+ ~*~ 
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® — 03 n=n 
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iK± 
e mn 


m 

r~ e + k± e 

p s $ mn (p 


28 


mn 


-im0 ~iK~ 2 
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/ 


dtu'6( u '- 


- i_ 

2tt 


oo 

E E a± 
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^ (F m p)‘ 1( p* + *3*0 


m 
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20 mn ( M m n P s) e ~ 1 (®^ + K+ ; 

v mn 
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APPENDIX D 


The geometry is shown in figure 1 6 , assuming for this discussion a 
stationary row. The blades are unstaggered, flat-plate, constant 
chord airfoils, so that e, = 1 and e =0. Letting z. . =0, then 

cp Z n • O t 

when points within the blade row are to be considered, the pressure 
spectral density resulting from the airfoils being dipole surfaces is 


P (r,u) 


CO GO 




im<{> 


m =-oo n=o j=o 




1K + 
mn 


(Z-C’) 


+e(€’ -Z’) 




I 


d€’ 


(Dl) 


with 



2 


(D2) 


and 


<± = - K± 

mn 2 mn 


(D3) 


The ^-component of the velocity spectral density associated with 
this pressure field is, from equation (A18) . 




V (p,<i>,Z,u)) 
<p 


1 ^ i — (Z-Z') ~ 

v, (p,<j>,Z=-<»,u)) — - / e ^ 1 9 P (p»$»Z* , 0 )) 

* V p 9 $ dZ 

—00 


2 00 00 N-l 2 


,u) + jjg ifijjp E E Er«.(w) 

ra^-co n=o j=o x ' 

im («- ^f]f - 

^ £ mnj ( « >“> 4 ra , ( «'- Z ) «' (D4) 


with 


f . 
mnj 





(D5) 


and 


S^CS'.Z) -«D e 1 "''' 5 ’ + 0 O'- 2 ') ^ "'’'j 


dZ* 


(D6) 


The boundary conditions at the airfoils are that this induced velocity 
cancels the "incident" velocity normal to the airfoil, 


u , FOR k = 0,1, 


N-l, 


_ 2 irk 

*k -t , 


u^(p ,Z j'O)) = u (p,4> k ,Z,w) 


v * C p * ^k ’ z ’ 


(D7) 
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with 


c < < c 

n < p < 1 AND - ~ Z — y 


and that the pressure difference or dipole density vanishes at the trail- 
ing edge (the Kutta condition) . These conditions turn equation (D4) into 
the system of equations for f 

mnj 


CO OO 2 

m 


-u^p.z.u,) = v +k (p, ?.-»,») + E 


<R 

iMp ^ $ m 

m=-a> n=o mn 


OVnO 


N-l 


■/{ 2 

-1 1 j=0 


2irim 


k-j 

i m 


(D8) 


N f. (£',&>) > tU\Z) dC 


mn 


th 

for k = 0, 1, ... N-l, and p,z on the k airfoil. The state of affairs 
in solving this problem seems to be as follows: writing equation (Dl) 

in the form 


p (r,to) 


N-l 

E ffh 0W“) 


j=o 



(D9) 


with da a surface area element, and using T = T™ + T_ (see appendix B), 
a S a r r IS 

with T the free-field propagator and T the "duct wall effects" func- 
r r B 

tion, gives 


P 


P + P 
FF B 


(DIO) 
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with p™, and p„ given by equation (D9), with and r_ in place of T, 

r r d * r r d 

respectively. The term p is then neglected for the purpose of solving 
for F_. . Further, the three-dimensional problem remaining is converted 
to a two-dimensional cascade problem by taking the airfoils to be 
infinite and parallel (see Kaji and Okazaki, [ref. 45] and Mani [ref. 46]). 



APPENDIX E 


The phase factor 


-in N, — TAN tj' 
e lp 


(El) 

in equation (3.1.9) accounts for the relative phase of the wake at the 
downstream row for different span positions. If tan ip is not propor- 
tional to p, then the wake will be skewed, i.e., nonradial, and the 
relative phase will be different for different span positions. Three- 
dimensional flow effects might also cause skewness in the wakes. A 
crude model to account for some of this extra skewness is the following; 
assume that the major three-dimensional correction will be in the tip 
region and consider the linear adjustment 


iTAN t-^TANg, (E2) 

where <p and this linear phase adjustment is illustrated in figure IT. 

This correction is included in the computer subprogram but has no 
other Justification than that it was a part of a Btudy to understand the 
effects of wake skewness on tone noise. The correction was not itself 
developed from theory. 




APPENDIX F 


The equations for the mode amplitudes of the four potential flow field 
interactions are very similar. They only differ in signs and indexes. 
All four interactions are therefore represented by the same set of 
equations. Modifications to the signs in the equations are done with 
switches C3 to C14, defined in the table following the equations. 
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K 
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(FI) 
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M - M + 

I,K Z,K 
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(F7) 
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(F9) 
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mn e Z 


<R (v P) V 1 
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The following table gives the values of the switches required in the 
previous equations to make the equations represent a specific interaction. 
For an interaction specified by a column in the upper box, the corres- 
ponding column in the lower box gives the appropriate value of the 
switches. The abbreviations and indexes used in the table are as 
follows : 


Component 

Index 

Abbreviation 

Inlet stator 

1 

IS 

Rotor 

2 

R 

Outlet stator 

3 

OS 


Velocity-inducing 

component 

IS 

R 

OS 

R 

IS 

R 

os 

R 

Li f t-producing 
component 

R 

OS 

R 

IS 

R 

OS 

R 

IS 

Subscript of velocity- 
inducing component, 
K2 

1 

2 

3 

2 

1 

2 

3 

2 

Subscript of sound- 
producing component, 
K1 

2 

3 

2 

1 

2 

3 

2 

1 

Sign of k 



<o 



>o 



K 

Z 

a 

Z 

a 

Z 

a 

Z 

a 

C3 

-1 

1 

1 

-1 

1 

-1 

-1 

1 

C6 

-1 

1 

-1 

1 

-1 

1 

-1 

1 

C7 

-1 

-1 

-1 

-1 

1 

1 

1 

1 

C8 

1 

-1 

-1 

1 

-1 

1 

1 

-1 

C9 

1 

-1 

-1 

1 

1 

-1 

-1 

1 


Cll 


-1 


1 






APPENDIX G 


Figure 18 illustrates the rotor blades as a two-dimensional cascade and 

defines the blade coordinates. It is seen that a point in the velocity 

th 

field is given in terms of the blade-fixed coordinates of the j blade 
by 

Z = Z| COS y + Y! SIN y + Z„ _ 
j J M.C. 


Y = - Zj SINy + Yj COS y + jb + Ut - £2ttp 


(Gl) 


where the i. 2 irp term results from the periodicity on 2 rp of the distortion 

(this must be modified if there is preswirl to the rotor) . Assuming the 

two-dimensional representation of the distortion, then the magnitude of 

tli th 

the velocity disturbance at the j rotor blade on the i revolution is, 
at Yj =0 
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FIGURE 18. - CASCADE REPRESENTATION OF ROTOR ALONE 
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The sum over all the periods gives the velocity 


W 



(G3) 


From the Poisson sum rule, 


?). 


SL= -DO 


n= -oo 


(GA) 


so that, using b = 2irp/N and U = M^p, 
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Assuming the distortion is frozen convected with the mean axial 
velocity in the duct, then 


W(Y,Z,t) = W(Y,Z - M z t) 
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giving 
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and, from equation (G5), 
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Finally, from the periodicity in y, 
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so that 
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This gives for W, letting = 2iTj/N, 
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When the distortion is steady. 
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so that 
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with 


v = 
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nM, 

~V 


T c 
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and 


V = 


M p 
SINy 


(G17) 


(G18) 


If the distortion is strictly axial, then the spectral density of 
the normal component of the velocity perturbation on the j ^ blade 
required in equation (2.3.1) is 


A 

u . 
3 


SIN y 


ou 

/ 


W(Z”=0) 


+ iwt 
e 


dt 


SIN W n eln< *3 6(w +nM t ) 

n =-oo 


(G19) 


When the distortion is nonsteady, the assumption is that it is the 
quasi-steady disturbance resulting from the convection with the flow of 
stretched eddies. Thus, the factorization can be made 


W = W(4>)AU-5 0 ) 


(G20) 


where : 


5 " z - M z t 

? o = Z M.C. " V* 

x = time at which the center of eddy is coincident with the midchord 
plane of the rotor 



Choosing to represent A(£-£ ) by a Gaussian distribution, 

- (£ -Sft ) 2 

A (5 - C 0 ) = e 2L 2 


(G21) 


with L the "axial length scale" of the eddy, then 


A(A) = ^7" Le lXCo e 2 


2 2 
il r 


L — * oo 


2tt6(A) 


(G22) 


If the eddy velocity disturbance has both axial (W z ) and circumfer- 
ential (W.) components, then the spectral density of the normal component 

♦ th 

of the velocity perturbation on the j blade is 


Uj = J | SINy W z (Z'=0) + COSy (Z=0)| 


iiot , 
e dt 



(G23) 



From equations (G9) and (G10) it is seen that when the distortion 

is nonsteady, the harmonic gusts at the blades have phase velocities 

differing from the relative mean flow velocity, i.e., the harmonic 

gusts are not "frozen” in the fluid. It is only when A /2 ttM and 

z z 

A^/2irM z approach being delta functions that these gusts became frozen 
convected. 



APPENDIX H 


Figure 12 illustrates the oblique gust geometry. Filotas (ref. 36) 
developed an approximate solution for the lift response of a thin air- 
foil to such a gust. His results were in a form containing an infinite 
series , 


T(h,</0 


1 

J h + F(h,<M 


I e (h 2 ) + I 1 (h 2 ) 

J ° ( h r ih D +iJ i C h r ih 2) 


(HI) 


where h^ = h sin ip and h 2 = h cos \p (equation [64] of ref. 36 with his 
k -*■ h and 0 -*■ ij») . In this equation, 1^ and 1^ are modified Bessel 
functions of the first kind, and are Bessel functions, and F is 
given in terms of an infinite series involving modified Bessel functions 
of the second kind (equations [21] and [22] of ref. 36). A generally 
good approximate form to equation (HI) was also developed: 


T = 



yj l+irh (l+SIN 2 iJ> + ?rh COSijO 


(H2) 


Both equations (HI) and (H2) have been programmed for use in sub- 
program BBCAA. Figure 19 is a plot of |t| versus h for different ip 
values. When <p = ir/2, the gust is no longer oblique and T then 
approximates the Sears function. 












APPENDIX I 


This appendix contains the list of equations used in the analysis that 
are referenced in the FORTRAN variable dictionary and the subprogram 
descriptions . 


1. 



2 2 
U-M) / 

mn 


2 . 


I T 


■ (o W MtS „n 


mn 


1-M 


, ON M 

3 * y ^ — - = e 

mn f ~ B 


Vl-M 2 


4. J* (nx) 

jiw- v? /— y :( x) = o 


m- ' Y' (nX) m 
m 


5. 


( ( n x ) > 

Y > x) { j ; <x) - ^ (x) } 


= 0 


J' n) 

R (y m p ) = J (pp) - - v ? , mn y (y P ) 

mn m mn Y (y q ) m ' mn 

m mn 


Nmn " [ 2 ^ ” „ 2 ) " 2 " „ 2 ) (v 

mn mn 


1 

2 
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8 . 


<R 

m 


(up) 

mn 


R (y p) 
m mn 


N 


mn 


9. 


A~ 

mna 


( CONSTANT | 

n sub 

1 f 

| FACTOR j 

l\S i 


bj 

AVERAGE OF 
NON-OSCILLATORY 
FACTOR 

3 a l 



OSCILLATORY 

FACTOR 


10 . 


V p > 




'Mj(p) - M 2 (p) 


+ Vm 2 ( P ) - M 2 (p) 


y 


i 

4 


+ Mj(p) 


11 . 


y(p) 


ARC COS 


M z (p) 

w 



13. 


e Z (P) = 




(P) 


14. 


1 mna 


(P) = 


_ C(p) 


+ 

K - • 

mn 


e,(p) - 


m 


e z (p) 


'] 


15. S(v) = 


^ ^Ho 2) (v) -iH^ 2) ( v) 


) 


134 



16. 


V*> ~ J o(v) + iJ 1 (v) 


17. 


r(2) 


T(v) = 


-H 


(2) 


(v) + 1 h| 2) ( v) 


(v) + iH^ 2 ^ (v) 


18. 


J (X) = J 0 (X) + i j 


(X) 


19. 


F (X) = T (X) • j 


[ J 00 ' i V x) ] - [. 


J(x) “ x J 1 (X) 


20. 


F f (v) - F(v) + A j( v) 


21 . 


L<V) = ^ ’ S(V) + b 2 *V v > + b 3 F f (v) 


22 . 


L ’(v) - b 2 . S(v) 


' J W + b2 + 


23. 


K«) * f J (^.)- »c> 
1 v ' 


MAX 


(~l) j . 


j-i 


Jj(v) 


f J j+l ( Kma ) +J j-1 (* 


inner) 
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24 . 


H^ 2 ) (v) 


Kjv.A) 


( . 1 ^ 

I [jo(A) - i J^A)] ' H^ 2 ) (v) + illi 2) (v) 


+ iv/A J 1 (A) 


H^ 2) (v) 


H^ 2) (v) + iHj 2) (v) 


Jo (A) -i J^A)] 

v 1 




IF vA =j= 0 
IF v=f 0, A = 0 


IF v= 0, A+ 0 
IF v= 0, A= 0 


25. 


T(v,0) = 


r tt 0 ( i+ coso ) ”1 

-iv SIN0 -- A -j 

L 1+2-rrv / 1+ j COS0) J 


'l+TTV 


1+SIN 0+rrv COSO 


)]* 


Additional equations for the viscous wakes interaction model: 


26. 


COS(ip) = 


M 2Z (P) 

Mie(p) 


27 . 


6(p) = y(p) + <Kp) 
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28. SIGOL *Y 0 • N 

T 1 “ 2 TipCOS *p 
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Average of 


non osc: 


factor 


SIN 6 

M A • A 

E,1 p-COSi|> a 


COMPACT 

SOURCE 


• S(v c ) - a C0T6 *F a (v a ) - f* COTg- F f (v 


•■] { 


me <j> + ' 

— — + K e 
p mn ^ 


Average of 


non osc 


factor 


illatoryU V C 2 {t^)' V 


™ SIN g . . 

M E,I p-COSi^ * A °’ A o 


i NON-COMPACT 
SOURCE 


— i. + Kiel • ["s(v )• J (kT \ - COTg -a- J fv + COTg fF'(v_) 

|_p mnz o \ mna J \ o mno I f o 


36 . j oscillatory J , .Wy _ SIGOL . Bl 


factor 


(“m" 0 ) 


Additional equations for the potential flow field interaction model: 


F K2 (p) = TlC K2 


^ M M,K2^ P ^ ^ A K2 ( - P ' ) + ^2^^ 


i th 

where A^ = Glauert coefficient of j order of component 
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a K,Kl (p) ' 


C 6 C K 1 (P) N K 2 - C ? K 


.n +1 .n -1 x 

A K 2 (p) ' A K 2 (p) 
8 n,K 2 P IA° 2 (p) + A K 2 ^ P ^ 


where 


A K 2 (p) ’ A K 2 (p) ’‘ **• A K 2 (p) ARE INPUTl 


A K2 1(p) = A K2 2(p) = 0 AND n =l,2,“*,N+ 1 


h - C 12 KN K2 C K2 (P) “13 [l ' C 14 e K2 <P) J 

K 2 ^ P; 2 p e 


H k k 2 (p) - J. (\ 2 < p >) + £ <5ll-f)X,K2< p )- J „C\2< P >) 


r, _ V.jlp) 

d K ,KI <p) ’ C 8 S K1 (P > + C 9 KS K2 [C b '' p ) TAN0 Kl( p ) + ' 2^ ^(p) 


43, „ , , r 

‘• K1 ' 6 2 M M,K 1 <P) 
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44. 


K e 1C 3 [I - 6 K1 (0) ] 


X k,K 1 <p) - c 6 


2p 


45. 


Constant 

factor 


N. 


K1 


48 


-iK“ Z _ „ 

■ e mn Kl*2ir 


mn 


46 . 

Non-oscillatory 
factor 


47. 


Oscillatory 

factor 


" ^Sl,Kl/ P ^ P K2^ P ^ 


a K,Kl (p) H ,c,K2 (p) 


e- ld «,Kl (p) 



Additional equations for the rotor-steady velocity distortion inter- 
action model: 


48. 


2tt 

W £ ^ = J W (P » < t > ) e i ^ < * > d< t> 


where 


^ 1-Ap COS 4>- (ApC0S4.-l) 2 - (a 2 -5(p 2 -i) 


w(p,40 = 



49. 


W 


,( p ) = \ a-^P) * 1 q 


50. 


W £ (p) = 


\ (a £ (p) + ib £ (p)) IF £>0 


\ (a (p) - ib z (p>) IF £<Q 


51. 


C (p)* M 

V p) = £ 2M m ( P ) 


52 . 


COT ( (3) = 


M z (p) 


& 


(P) - M z (p) 


53 . 


„ -IK* Z„ 
„ i N e mn R 
Constant / R 


factor 


43 


mn 


54. 


Average of 

non-oscillatory 

factor 


= C r 


(^r} • M z [t 4 + K± - e - 


R \da / ^ "Z 

COMPACT 
SOURCE 


p mn Z 




) - d -COTS-FaCVj) - f C0T3 -F f (v 4 ) 
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55 . 


[ S(v ) • J [k + 1 - COTS • a -J L + k 1 \+ COTS* f- Fi 
V \ mna J 1 £ mno I f 


f (v £ 


Average of 

1 / dC T \ [ me <U + 

non-oscillatory 

= C ( -j — — — )• li. • M • SINg- — ^ + K~ e 
R\da / >1 Z p mn Z 

factor 

NON COMPACT L J 


SOURCE 


56. 


\ Oscillatory I _ w 


(P)‘ « 


I factor 


m 




Additional equations for the rotor nonsteady velocity distortion 
interaction model: 


Average of 

non-oscillatory 

factor 



+ K~ e„ 
mn Zl 


58. 


1 Constant 
j factor 


N_ 

4 


-iKT Z 
e mn 

^mn 


\ Oscillatory 
j factor 


FACTIN4 


60. 


T. 

3 
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' 1 / X \ 2 


61. 


P. = U. 
3 J 


2 

e 

IV 

IF 

BT. > 10 
J 



IF 

.1 < BT . < 
J- 

2BT. 

3 

SIN(Bt) 



^2 IT * 

Bt 

IF 

BT. < . 

3 

T. 

e = — i 
3 

} _<& 
= J * 2 

-B 

) 2 

* COS (lOT ) da> 



If ARMISC(25) = 3: 

62. FACTIN4 = <R 


m 


(P mn P)* Mjj(p) -SIN(y) (p ) -F 1 (p) - S 2 £ (P) *F 2 (p)J 


C 2 (P)* 

'A'"' ~ 2~%&T 


' o \ K / * 

63 * v.(p) =£* 


64. 


[ F i (p) ] 


- S(v £ ) - a(p)- COT(y)- F a (v £ ) - f (p) * COT(y) • F f (v £ ) 


, COMPACT 
SOURCE 


65. 


[ F i(p>] 


NON-COMPACT 

SOURCE 

+ 


S( V’ J Gmo) " “ WC0I <’)' J CV K m IW > £(o) ’ CO ' r( ' f) ' F f ( 'V 


66 . 


i 


f 2 (p) 


= COT(y) *S(v £ ) + a(p)-F a (v £ ) + f(p)- F f (v £ ) 


COMPACT 

SOURCE 
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67. k(p) -COTtlf) -SCv,) .J(J^ + d<p).J(v t +K^ o ) 

L J NON-COMPACT 

SOURCE + f(p)- F £ ( V £> 


68 . 


g. (P) 


(P-R) 2 
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A = , KMAX = - 7 ^ +1 

K A 


h(p,k) =^Vk 2 + ± SIN(y)' 


72. p ( P ,k) = TAN — p SIN (y) 
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APPENDIX J 


The zeros of A (y) are symmetric under a change in sign of m, i.e., 
m 

for 


A 

m 


(y) = J'(y) Y' 
m m 


(yn) - J' (yn) Y' (y) 
m m 


(Jl) 


A -m (y) = A m (lJ > 


(J2) 


Hence, the zeros are computed for A j | (y ) . All the zeros are nonnegative. 
They are computed by first approximating the zeros for m = 0 as follows: 


1) For q = 0, use formula (9.5.13) of reference 30 

2) For .2 = q < 1, use formulas (9.5.28) and (9.5.31) of reference 30 

3) For 0 < q < .2, use quadratic interpolation with the values for q = 0 
in item 1 above and q = .2 and .3 in item 2 


Using the approximation formulas in item 2 restricted to q = .2 
yeilds a one-to-one correspondence with the eigenvalues (zeros) where 
the first several values are poor approximations. The approximations 
in item 2 get worse as q decreases and for q < .2 no longer give a one- 
to-one correspondence. For this reason, interpolation procedure 3 is 
used to obtain better values. For more accurate values for the zeros, 
an iteration method is applied with above approximations as starting 
values. Reference 47 gives an iterative procedure for computing a real 
zero of a real valued nonlinear function f(x) in which a rational 
function is fitted through previously computed values. The interaction 
formula is 


n+1 


X n + 


(x -X .Vx -x n f (f ,-f 

V- n n-l.yV n n-2y n V. n-1 n-2 J 

rx -X ff 0 -f > f . + (X -X ( f -f .M 9 

V. n n-ly v n ~2 nj n-1 V n n-2 J V n n-1 J n-2 


(J3) 


146 



where f if f(x ). Here, f is Ai i and x is u. One of the three 

n n [ m[ 

starting values required is chosen form the approximations discussed 
above, with the other two that approximation ± s where e is chosen to 
"surround" the zero. 

Applying the iteration until f(x^) is small provides zeros with accuracy 
limited by the evaluation of f itself. In this case this accuracy is 
that of the Bessel and Neumann function evaluators. Do this for m = 0, 
use the m = 0 zeros as starting values for the m = 1 zeros, etc. 

Examples of the zeros used as eigenvalues in the analysis are plotted 
in figures 20, 21, and 22, each figure corresponding to a different 
value of n • 
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ZEROS OF A m ( M ), F = Fmn> n ~ 0, 1, 2, . . 
NON-DIMENSIONAL INNER DUCT RADIUS V 
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